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ABSTRACT 



a ( 

This thesis concerns developing a program that takes an 
aerial photograph, and a set of Digital Terrain Elevation 
Data (DTED) that is defined over the area of the photograph, 
and generates a synthesized view that represents what a 
camera would see from a different location. The elevation 
data points are grouped into triangular panels that are 
projected to the reference image by three dimensional 
transformation equations. Shading for the synthesized image 
is determined from the reference image. The pixels of the 
reference image that fall within a triangular panel are 
collected and averaged. When a new observer location is 
selected, the panels are projected to the new synthesized 
image plane. A z-buffer approach and a polygon fill 
algorithm were used to remove hidden surfaces of the 
synthesized view. 

This program is tested on both artificial and real data. 
Other characteristics and performance measurements of the 
program are also analyzed here. The quality of the syn- 
thesized image from real data was affected by the low 
resolution of the terrain elevation data, and yielded less 
desirable results than could be expected of a higher 
resolution terrain model. 
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I . INTRODUCTION 



A. COriPUTER I NAGE GENERATION FROM AERIAL PHOTOGRAPHY 
The main objective of this study uas to develop a 
program that takes a digital photographic image and a File 
of terrain elevation points defined over that image as 
input, then produces as an output a synthesized perspective 
view. The synthesized view is a rotated 3-dimensional (3D) 
perspective representation of the original photographic 
image. The main application of this study is to generate a 
different perspective of a terrain model. This may be used 
to generate different views that a pilot of an aircraft 
could expect following different Flight paths through the 
same area. Further study may make it Feasible to generate 
synthesized images fast enough to simulate a real time image 
display of a flight for a mission briefing or to be used as 
a training aid. Another application could be for training 
men on optically guided missiles. With high resolution 
images, a flight path through a battlefield could be 
simulated that would have all the visual characteristics of 
an actual flight without the expenditure of a live missile. 
The generation of a shaded image as a 3D picture provides 
unique problems for 3D graphic displays. The data which 
comprises a photographic image consists of an array of 
pixels, each of which has a defined grey level or shade. 
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There are 256 different levels of grey that may be assigned 
to a pixel . This study concerns taking a photographic 
perspective image and a 2-dimensional (2D) array of eleva- 
tions defined in a grid covering the area of the image as 
inputs. The grid of elevations, called the terrain model, 
is geometrically related to the photographic image through a 
perspective projection transformation that equates the world 
coordinates of the elevation points to the object coor- 
dinates of the image. A synthesized image from a different 
observer location is then generated. The new synthesized 
view should approximate what would be seen by a camera from 
the new observer position. 

The differences between the original and synthesized 
images will be affected by the resolution of both the 
photographic image and the terrain models. Higher resolu- 
tion of the original models will result in a closer approx- 
imation in the synthesized image. Another complication or 
ambiguity arises when details which should show up in the 
synthesized view were not present in the original image. A 
method must be devised so that it can fill in areas which 
become visible in the synthesized image that were hidden in 
the original reference image. The solution to this hidden 
surface problem is further addressed in the discussion of 
the grey scale referencing a Igor 1 thm . C Ref . 13 
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B. TERRAIN ELEUAT I ON AND PHOTOGRAPH AS INPUTS 



In this study a high altitude aerial image of Moffett 
Field, California was used as the original reference photo- 
graph. The photographic image was supplied by the Defense 
Mapping Agency (DMA) and had a resolution of approximately 1 
meter per pixel . The terrain model corresponding to the 
reference image was provided as Digital Terrain Elevation 
Data CDTED) by the DMA and consisted of elevation points 
taken every second of a degree change in latitude and 
longitude. This gives an approximate resolution of 30 
meters per elevation point in a north-south direction and 23 
meters per elevation point in an east-west direction. 

The synthesized view was restricted to a northerly 
direction which simulates an aircraft flying from south to 
north with the image plane perpendicular to the direction of 
flight. To allow for different flight patterns would 
require development of an algorithm that would provide for 
rotation of the image coordinates which is beyond the scape 
of this study . The main idea is to generate a synthesized 
view that is rotated from the original photograph by 
approximately 90* and explore the concepts of the algorithms 
required to da this. Although speed was not a major issue, 
the size of the terrain model was limited to a 50 x 50 grid 
array, or 2500 data paints, to minimize the time far 
synthesized image generation. 
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C. ALGOR I THM ISSUES 



1 . Greu Scale Referencing 

To determine the grey scale values of the pixels 
that make up our synthesized view, the terrain model is 
first divided into triangular panels. The vertices of the 
triangular panels are then mapped into the original OHA 
image using the perspective projection transformation that 
projects georectangular coordinates into reference image 
coordinates. The pixel grey scale values that fall within 
the projected triangular panel are then averaged. The 
average grey scale value is permanently assigned to that 
particular panel . When the synthesized image is constructed 
the triangular panel is mapped to the new image view and 
filled with the assigned average grey scale value. In this 
way the sample of pixels that fall within the triangular 
panel are mapped from the original reference image to the 
synthesized image. This method of mapping the triangular 
panels to the synthesized image also salves the problem of 
assigning a grey scale value to hidden surfaces of the 
reference image because they are automatically assigned the 
value of the surrounding pixels. Since the average grey 
scale value for a triangular panel is dependant upon the 
resolution of the terrain model, the grey scale value 
assigned to the hidden surface will also be affected 
CRef. 13. 
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The smaller the triangular panels are, the smaller 
the area that must be collected and averaged in the original 
image. This means a much better synthesized vieui can be 
constructed that contains more of the attributes of the 
original image. For this reason the resolution of the 
terrain and reference image model is very critical to 
obtaining an accurate synthesized view . Using the resolu- 
tion of the terrain and image models used in this study, the 
approximate number of pixels that must be averaged in the 
reference image for each triangular plane mould be 1/2 C30m 
x 23nOCl pixel/m} - 345 pixels. This is very coarse and 
does not allow for optimal generation of the synthesized 
view . 

2 . Hidden Surface Elimination 

There are surfaces that may be discernible in the 
reference image but become hidden in the synthesized view. 
The z-buffer algorithm was used to accomplish the hidden 
surface elimination. The z-buffer is an array that 
contains the depth or distance to the observer location for 
each pixel that is to be visible in the synthesized image. 

As each triangular plane is mapped to the synthesized view 
the location and depth of each pixel within the plane is 
determined. The depth of the pixel to be written at a 
certain location is compared to the depth of any pixel that 
may have been previously written to the same location. If 
the depth or distance of the new pixel to the observation 
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point is shorter than the previous pixel, its depth is 
written to the z-buffer and the grey scale value of the 
pixel is placed into the synthesized image. If the depth is 
larger, no updating occurs and a new pixel is obtained in 
the process. The z-buffer works very closely with the next 
algorithm to be considered. CRef. E, pp . E65-E673 
3 . Polygon Fill Algorithm 

Screen coordinates are generated for the three 
vertices of each triangular panel as it is mapped to the 
synthesized image. Screen coordinates are designated as IP 
and JP values with the IP values representing the columns 
and the JP values the rows. The location, IP,JPC0,03, 
designates the upper left hand corner of the screen and the 
maximum screen coordinate, IP , JPC51E , 51E) , the lower right 
hand corner. Pn active edge list CPELD is generated by 
computing a line between each of the translated vertices. 

For each line, the IP coordinate corresponding to the 
maximum JP value of the line, the amount IP changes for each 
one unit step of JP, and the total span of JP are stored 
into the PEL. The three lines generated for each translated 
triangular plane will form another closed triangle. By 
using the parameters stored in the PEL the location of 
pixels enclosed by the translated triangular plane can be 
determined, and the corresponding array points within a 
frame buffer are changed from a 0 to a 1 . 

CRef. E, pp. 75-793 
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After all of the enclosed pixels have been marked 
within the frame buffer, the buffer is scanned row by row. 

If a value of 1 is found, then the depth is calculated for 
that point and compared with the depth value stared in the 
z-buffer. If the depth value is smaller, that pixel is 
located closer to the observer location and the grey scale 
value for that pixel is written to the synthesized image 
file. As can be seen the fill and hidden surface algorithms 
work together to generate the new image. The implementation 
of these algorithms are explained in further detail later. 

In Chapter II there will be a discussion of basic 
photographic geometry to develop an understanding of the 
transformation equations necessary to map object coordinates 
into image coordinates and for image plane rotation. 

Chapter III will detail program considerations based on 
image and elevation data as well as the algorithms used to 
generate the synthesized view. Chapter IU will discuss 
passible ways to improve the transformation program and 
discussion of topics for possible further study. An outline 
of the program is contained in Appendix A that gives a short 
discussion of each subroutine as to its purpose, input and 
output, modules that called, and modules that reference the 
subroutine. Appendix B contains the entire 3D transforma- 
tion program . 
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I I . PHOTOGRAPHIC GEOMETRY 



A . BACKGROUND 

To understand many of the concepts used in this study, a 
basic background in photographic geometry is presented. The 
relationship between the image space and object space is the 
basis For many of the equations that help to generate the 
reference and synthesized images For visual display . The 
objective of this chapter is to present the concepts that 
allow the transformation of 3D objects to a ED image and the 
parameters evolved. 

1 • Perspective and Parallel Pro lection 

A parallel projection is a projection m which the 
projection lines from the object to the image plane never 
converge. When an object is viewed by parallel projection, 
its size would never change as the camera is moved closer or 
Further away. In contrast, a perspective projection has all 
the projection lines From an object converge to a perspec- 
tive center. A perspective projection imitates how we see 
things. An example would be a picture of railroad tracks. 
The tracks would appear to become closer together when 
Further away From the observation point. In a parallel 
projection the tracks would be the same distance apart along 
the entire length. [Ref. 3, pp . 133-134] 

Since a camera is generally designed to photograph a 
rather large area, it involves perspective projection. The 
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camera view represents what an observer would see standing 
at the same location, and the images generated are perspec- 
tive images. This means that the equations used to trans- 
form the object space into the image space must be perspec- 
tive transformations . 

2 • Image Coordinate Space 

The image plane is the plane of the photograph to 
which the object points are mapped. It has a 2D coordinate 
system to which each point of a 3D object is translated to 
CFig. 2.1). The indicated principal point CIPP) is the 
center of the image plane and has the coordinates of 
Cx,y,0). The x, y, and 2 axis represent a right handed 
plane and the perspective center CL) lies along a line 
parallel to the z-axis; then a perpendicular line is drawn 
from L to the image plane. The point at which this perpen- 
dicular line intersects the image plane is called the 
principal point Co) . This offset of the principal point 
from the IPP is compensated for in the transformation 
equations by xo and yo. The focal length of the plane is 
defined as the distance from the principal point to the 
perspective center. For the image plane coordinate system, 
each object point CA) is graphed to a corresponding image 
plane point Ca) located at (xa,ya,0), and the perspective 
center or focal point is located at Cxo,yo,f). CRef. 3, pp . 
135-1363 
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In generating a synthesized view, the perspective 
center may be placed at any location desired with reference 
to the image plane. It is therefore desirable to select a 
point along the z-axis such that xo and yo become 0. This 
will decrease the number of calculations required in 
generating the synthesized view. 

3 . Obiect Coordinate Space 

The observer location and each object point position 
is described in world coordinates, called georectangular 
coordinates, of X, Y, and Z. The center of the earth is 
given as CO, 0,0), the Z-axis points directly to true north, 
the X-axis points to the intersection of 0* Latitude and 0* 
Longitude, and the Y-axis the intersection of 0* Latitude 
and 90* E. Longitude. The principal or focal point that 
would describe the observer location in georectangular 
coordinates is XL, YL, and ZL . Each object point CA) 
located in the object space is identified by XA , YA, and ZA 
as shown in Figure 2.2. CRef. 3, p. 13SJ 

B. iriAGE PLANE ROTATION 

To align the x, y, z coordinates of the image plane to 
the desired viewing direction requires rotation about the X, 
Y, and Z axis of the georectangular coordinate system. In 
general to transform one 3D coordinate system requires a 
matrix multiplication of the form A = CHUB. The A repre- 
sents a vector in the image space with x, y, and z coor- 
dinates, and B is a vector in the georectangular coordinate 
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z 




Fig. 2.2 Object Plane Coordinates (Ref. 3, p. 136) 
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system with X, Y, and 2 components. 



This may be written as 



X 
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This maps the vector X, Y, 2 to the image space x, y, z. 

The n matrix is derived from the definition of the direction 
of a vector A given by 



= Case; 1 + CosB j + Cost k 

A 

Ax Ay Az 

A 3 A A 

= i + + j k C2.3D 

A| |A| | A 

A A A 

where i, j, and k are unit vectors of the particular 
coordinate system in which vector A is contained. The 
quantities <x , B, t are the angles that vector A makes with 
the x, y, and z-axis respectively. Since there are three 
vector components in the X, Y, 2 system, each one must make 
its own transformation into x, y, and z and thereby forming 
the 3x3 matrix of 11. To translate x, y, z into X, Y, 2 
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the inverse of the M matrix is taken giving B = Cnif'A. 

CRef. 3, p. 1392 

If we can define three orthogonal vectors in georect- 
angular coordinates as R, S, and T that would describe the 
desired viewing position of our image plane, the fl matrix is 
easily derived by 





Sx 


Sy 


Sz 










l S l 


l S l 


l S l 










Rx 


Ry 


Rz 








n - 


— 




— 




C2 .45 






I R I 


l R l 


I R I 










Tx 


Ty 


Tz 










_ |T| 


l T l 


|T I_ 








Generally the image plane rotation is expressed 


in omega 


Cw) , phi COD, and 


kappa C k D . 


The uj is the rotation about 


the X-axis, 


the 0 


is 


rotation 


about the Y-axis, 


and k is 


about the 2 


!-axis . 


If 


we rotated first about the 


X-axis by 



an angle of w radians, the terms of the M matrix would 
equate as follows 

CosxX = CosCO*) “ 1 

CosyY = CosCwl 

CosyZ “ CosC90*-w) = SinCw) 

CoszY = CosCw + 90*) ” -SinCwl 
CoszZ = CosCw} 

All other terms equate to cos 90* • 0, therefore the fl 
matrix becomes 
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1 
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0 



n = 



0 CosCuj 3 SinCw3 
0 -SinCu3 CosCuj3 
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Similarly a rotation 91 about the Y-axis produces 
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and for a rotation k about the 2-axis 
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By multiplying all three matrices together me derive the 
overall n transform in ui, 91, and k, 



n 



CosC 91 3Cos C k 3 Cos C u 3 Sin C k 3 + S in C u 3 SinC 0 3 Cos C k 3 
-Cos C913SmCk3 Cos C u 5 Cos C k 3 - SinCu3 SinC 91 3SinC k 3 
SinC913 - Sin(u))Cos(J]) 



SinCw3SinCk3 - CosCu3SinC 91 3CosCk3 
SinCw3CosCk3 + CosC u3SinC 91 3 SinC k3 

CosCu3CosC 91 3 

C2.B3 



This is the general form of the matrix that maps the 
georectangular coordinates to the image plane. CRef. 3, pp . 
597-600 J 
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The general farm af the transformation matrix is used to 
initially map the elevation terrain model into the original 
reference image. The w, 9 , and k were supplied with the 
original DMA photograph and represents the rotation of the 
reference image plane with respect to the georectangular 
coordinates at the time the picture was taken. When we 
generate the synthesized view the image plane must be 
rotated to the desired viewing angle of the observer 
(northerly direction in particular} . This means that a new 
w, 9, and k must be calculated, or by defining new image 
plane coordinate axes in terms of the georectangular 
coordinates, we can calculate the terms of the M matrix 
directly using the preceding equations. This is discussed 
further in the next chapter. 
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III. ALGORITHM CONS IDERATIONS 



In this chapter an in-depth analysis of the original 
image and terrain data is presented. This includes how the 
data was referenced, the size of the data files, how the 
data was used, and how the elevation and image data compared 
with one another. The process of translating from the 
object space coordinates to image space coordinates and then 
to screen coordinates is considered, and the equations used 
are given. 

The referencing, fill, and z-buffer algorithms also are 
discussed in detail. How the data generated from these 
algorithms is used and put together to produce the syn- 
thesized view will be presented. Any problems that were 
encountered and the eventual solutions will be discussed in 
the appropriate section to which they pertain. 

A. REFERENCE IMAGE DATA 

The picture of Moffett Field supplied by the Defense 
Mapping Agency CDMA), was 4339 by 4997 pixels in size and 
came with both a left and right image. Only the left image 
was used to generate the perspective views in this study. 
Each individual pixel within the original image is desig- 
nated by coordinates I and J, where I is the pixel column 
and J is the scan row. The geographic northeast corner has 
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the I, J coordinates of C0,0), and the southwest corner 
C 4997 , 4999 ) . 

The image display devices used were capable of dis- 
playing images that were only 512 by 512 pixels in size, 
therefore, the original image was divided into appropriate 
blocks suitable For viewing called Frames. Each Frame 
contains an image that is 512 by 512 pixels. The coor- 
dinates of each Frame has a Four digit I_.Frame and 
J_Frame value, and is Further identified as a left or right 
CL or R) image. The I _Frame and J._Frame coordinates are 
designated in multiples of 512 which define the column and 
row location of each Frame. There are 10 Frame columns and 
10 Frame rows with assigned coordinates of 0000 through 
450B . To identify a particular Frame of the Drift image one 
first designates whether it is a Left or Right image and 
then give the I_Frame and J,_Frame coordinates. fts an 
example L05121024 would designate a Left image From the 
second column and third row. The first Frame L00000000 
starts in the southeast corner and the last Frame L45O046O0 
is in the northwest corner. The disparity between the 
starting location of the frame coordinates and the I and J 
coordinates of the original image must be compensated For in 
the equations that are used to determine the location of 
individual pixels within a Frame image as shown later. 
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1 . Object to Reference Image Transformation 

Every object paint is converted From its 30 X, Y, 
and Z gecrectangular coordinates into the ED x and y image 
coordinates using the A = CfIJB equation discussed earlier. 
The vector From the perspective center to each object paint 
is deFined by CXA-XL), CYA-YL) , and CZA-ZL) . This vector is 
mapped into the reFerence image plane coordinates oF x, y, 
and z. Since every vector in the image plane is directed 
From the perspective center or Focal point to the image 
point Ca), the z coordinate value is constant and equal to 
the negative oF the Focal point C-F) oF the camera. Using 
these parameters the equation becomes 
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a b c 
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uhere K is a scale Factor. From this transFormation the 
Following equations are obtained 

x-xo ■ K [a CXA-XL) + b CYA-YL) + c CZA-ZL 3 3 C3.2) 

y-yo - K [d CXA-XL) + e CYA-YL) + F CZA-ZL)] C3.3) 

-F * K [g CXA-XL) + h CYA-YL) + i CZA-ZL)] C3.4) 

Dividing the last equation into the First two and rear- 
ranging yields 
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X = xo 
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a CXA-XL) + b CYA-YL) + c CZA-ZL) 



g CXA-XL) + h CYA-YL) + i C ZA-ZL ) 



C3 .53 



y - y° 
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d CXA-XL3 + e CYA-YL) + c CZA-ZL) 



g CXA-XL) + h CYA-YL) i CZA-ZLD 
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uihere x and y are the ED image plane coordinates . C Ref . 3, 
pp. 141-1423 

The original parameters of the fl matrix mere 
calculated From the w, ID, and k that were given by the DHA 
with the original photograph. These represent the physical 
position of the image plane in relation to the georect- 
angular coordinate system at the time the picture was taken. 
The Focal point was also supplied and is particular to the 
camera that was used to take the original photograph. When 
the synthesized image is generated, the image plane is 
oriented to a position for the desired viewing angle, which 
means that the parameters of the fl matrix will change and 
must be recalculated. The steps used to determine the 
desired orientation of the image plane coordinates will be 
discussed later in this chapter. 

2 . Reference Image to Screen Coordinate TransFormation 

Once the x and y image coordinates have been 
calculated they are translated to I and J values of the 
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original image. This is accomplished using the affine 
transform which represents a ED into ED coordinate transfor- 



mation, The equation that accomplishes this is derived from 
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where Cl and CE are the values that translate the image 
plane origin to the I and J coordinate system origin. They 
are calculated by setting I and J to 0 and solving for x and 
y. To get the desired transformation of the image coor- 
dinates to I and J coordinates the inverse transform is 
taken. CRef.3, p.5931 
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Again the original J, k, 1, m, Cl, and CE were 
supplied by the DMA with the original image. Equation 3.B 
represents any general ED into ED transformation, and it is 
used both to translate the original image plane coordinates 
into the I and J coordinates of the reference image and to 
transform the image plane coordinates of the synthesized 
view into screen coordinates of IA, and JA. 

The screen coordinates were assigned the parameters 
IA, which represents the columns, and JA, which represents 
the rows. The point IA,JAC1,1) is mapped to the upper left 
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corner af the screen and IA, J A ( 5 1 2 , 512) is the lower right 
corner . 

The screen coordinates represent the location of a 
pixel within a frame. To convert I and J original coor- 
dinates to IA and JA screen coordinates one needs to know 
the particular frame one is working in and compensate for 
the difference in the starting location of the frame 
coordinates and the original image coordinates. This is 
accomplished by using the following equations, 

I A - Cl - I Frame) C3.9) 

JA - C 4999 - JJrame - J) C3.10) 

The I_Frame and J_Frame values must therefore be 
given to determine the screen coordinates within a desired 
frame. To allow flexibility in determining which frame 
image would be used to extract the reference image, an 
interactive input of the I_Frame and J_Frame coordinates was 
appropriate . 

3 . Sunthesized Image Plane Rotation 

For the synthesized views that were generated, the 
affine transform parameters Cl, C2, J, k, 1, and m that 
would map the newly rotated image plane into the screen 
coordinate system were selected. Since the image planes in 
the synthesized views were oriented for an observer looking 
north, the image plane coordinates were generated by cal- 
culating the z-axis, which points south, from two georec- 
tangular coordinate vectors calculated from two different 
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terrain data points along the same longitude line. By 
taking the difference between the X, Y, and 2 coordinates a 
third vector was formed that defined the image plane z-axis. 
The image plane y-axis was calculated by using only one of 
the terrain data points used to calculate the z-axis. By 
taking the negative of the georectangular coordinates of the 
terrain data point, a vector is produced that points 
downward through the center of the earth. This was the image 
plane y-axis. Once the z and y axes are calculated, the 
cross product of y cross z was used to calculate the x-axis. 
Figure 3.1 demonstrates the resulting image plane. 



N 




z 



Fig. 3.1 Synthesized Image Plane Coordinates 

This image plane coordinate system, from the 
observers perspective at C0,0,-f), would have the x-axis 
pointing left, the y-axis pointing down, and the z-axis 
pointing directly toward the observer. The screen 
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coordinates have the IA axis pointing right and the JA axis 



painting down . Therefore, the new values of Cl and C2 with 
IA, and JA equal to 0 were as follows: 

Cl ” naximum assigned x-image value C3.ll) 

C2 - 0 C3.12) 

which aligned the image plane x,yCo,o) to the screen 
coordinates of IA,JAC0,0). The J, k, 1, and m values of the 
transformation matrix were selected to scale the image plane 
to the screen . 

Having determined the three synthesized image plane 
coordinate vectors in terms of georectangular vectors, it is 
relatively easy to generate the M matrix parameters using 
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were S, R, and T are the georectangular coordinate vectors 
of the x, y, and z axes of the rotated image plane. This 
defines the new transformation matrix that will be used to 
generate the synthesized view by mapping the georectangular 
vectors of the terrain data into the new image plane. 
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B. TERRAIN DATA 



The Digital Terrain Elevation Data CDTED) supplied by 
the DP1A came as a rectangular grid of elevation data paints. 
Each elevation data point was recorded as an integer value 
in meters above sea level. If a particular elevation point 
was unknown, it was assigned a value of -32767 to assure 
that it would not be confused with any valid elevation data 
points. The rectangular terrain grid listed elevation 
paints every one second of a degree change in latitude and 
longitude. The southwest corner of the terrain grid was 
defined as being located at 37* 22’ 47” N. latitude and 
-122* 05’ 03” W . longitude. From this reference point the 
elevation data was laid out in 210 raws by 239 columns. The 
rows represented lines of constant latitude and the columns 
lines of constant longitude. With this information the 
northeast corner of the terrain grid was calculated as being 
located at 37* 26’ 17” N. latitude and -122 01 ’04” W. 
longitude . 

1 . Data Uerification 

The first problem was to compare how accurately the 
elevation data matched up with the original image data. 

This required taking specific elevation points that were 
known to match specific image points, then translating the 
georectangler coordinates of those elevation points to the I 
and J coordinates for comparison to the original image. The 
w, ffl , and k for the M matrix to transform gearectangular to 
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image plane coordinates and the parameters For the affine 
transform from image plane to original image I and J 
coordinates were supplied by the DMA with the original image 
data . 

To select specific elevation points for comparison 
required finding a method of distinguishing unique elevation 
patterns that could match specific objects in the image. 

The technique used was to visually display the elevation 
data as an image. The elevation image file was produced by 
assigning grey scale values to each elevation data point 
with the lower elevations receiving the darker shades and 
the higher elevations the lighter shades. When the eleva- 
tion image was displayed as shown in Figure 3.2, a distinct 
highway pattern emerged from which three intersections could 
reasonably be distinguished. The three intersection points 
Cshown as a, b, and c in Fig. 3.2) correspond to elevated 
roads that crossed one another in the original image. 

□nee the elevation points were selected for com- 
parison, the approximate row and column of the elevation 
data corresponding to the center of the intersections was 
determined. From the reference point of the terrain grid 
there is a 1 second change in latitude and longitude for 
each row and column which allowed the calculation of the 
latitude and longitude for each of the three reference 
elevation paints. The latitude and longitude for each point 
was converted to georectangular coordinates using a 
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conversion program, then transformed to I and J coordinates 
using the provided DPI A parameters as explained earlier. 

Each of the three elevation points mapped to within 10 
pixels of the original image in both the I and J coordinate 
directions. This equates to less than 1 second error in 
latitude and longitude which was deemed precise enough to 
establish the correlation between the terrain and image 
data . 




Fig. 3.2 Elevation Image 
2 . Elevation Line Drawing 

By selecting a smaller area of the DTED data, a more 
distinct picture could be studied. An intersection used to 
verify image and elevation correlation was selected to be 
used as the reference image. A smaller rectangular set of 
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terrain data that mould map to the intersection, plus a 
small section of the surrounding area, mas extracted from 
the reference terrain grid. This smaller set of terrain 
data mas taken from rams 71 through 79 and columns 172 
through 183 of the terrain model . 

To verify that this set of elevation points mould 
appear like the desired reference image, a line draming 
illustrated in Figure 3.3 mas created using a commercial 
graphics program called MDUIE.BYU CRef. 43. This program 
can generate a connected line draming from a set of eleva- 
tion data and alloms rotation as mell as magnification of 
the draming. To assure a sharp visual contrast, the 
elevation data mas magnified by a factor of 10 and the 
draming rotated to a useful vieming angle. 

The results mere not as desired mhich gave an early 
indication that the resolution of the terrain data may not 
be adequate enough to generate a synthesized viem that mould 
be a close approximation of the reference image. This mas 
further verified when the synthesized viem mas produced at a 
later time. Fin artificial set of image and elevation data 
mas used to verify that the transformation program func- 
tioned as desired before generating synthesized viems of the 
original reference image. The artificial elevation points 
mapped into the artificial image exactly may as the original 
elevation data mould map to the original reference image. 

The artificial reference image, shomn in Figure 3.4, 
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Fig. 3.3 Elevation Line Drawing 




Fig. 3.4 Artificial Reference Image 
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was of a small square structure resting on top cf a much 
larger square structure. Selecting a large abject Far an 
artiFicial reFerence image would minimize the eFFects oF the 
law resolution aF the elevation data. This produced more 
reasonable synthesized images that demonstrated the perspec- 
tive transFormat ion more accurately. The results oF the 
transFormat ion will be discussed Further on in this chapter 
aFter considering some oF the algorithms used to generate 
the synthesized view. 

C. SPECIFIC ALGORITHnS 

1 . Image ReFerencing Algorithm 

Due to the resolution mismatch between the reFerence 
image and the terrain data, a smaller subset aF the terrain 
model was selected. This would allow the transFormat ion oF 
smaller and more distinct images that could show the eFFects 
oF the program better. The desired terrain elevation points 
are extracted From the larger reFerence terrain grid by 
interactively selecting the proper rows and columns that 
deFine those particular points. The program allows up to a 
SO by 50 array oF elevation points to be extracted. This 
size was chosen to limit the time required to generate a 
synthesized view. From the rows and columns, the latitude 
and longitude is tabulated For each point, and is then 
converted to georectangul ar coordinates. The georectangular 
coordinates are written to a File in row order From leFt to 
right and From bottom to top. The First set oF coordinates 
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match the southwest elevation point, and the last set of 
coordinates the northeast elevation paint. 

The georectangular coordinates of each terrain data 
point is translated to I , and J original image coordinates 
using the camera parameters supplied by the Drift . They are 
further processed to Ift and Jft screen coordinates using the 
I_ Frame and J__Frame of the desired reference image. The Ift 
and Jft coordinates derived from the elevation points will 
now correspond to the Ift and Jft coordinates of the reference 
image. ftny four adjacent elevation points will define a 
rectangle which is divided into two triangular panels by 
defining a line connecting two opposite corners. Once the 
triangular panels of the elevation points are defined in 
terms of Ift and Jft coordinates, the corresponding pixel grey 
scale values of the reference image that fall within the 
same screen coordinates of the triangular panel are col- 
lected and averaged. fts seen in the example of Figure 3.5, 
the calculated average grey scale value is placed into a 
file along with the three specific elevation points that 
make up the triangular panel . This procedure is repeated 
until the entire terrain grid has been processed. Each 
triangular panel represents a sample or average grey scale 
value of the reference image. 

When the latitude, longitude, and elevation of the 
new observation point is input into the program, a new XL, 
VL, and ZL is calculated that corresponds to the new 
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perspective center. As explained earlier, the image plane 
is rotated to the desired viewing angle, which changes the n 
matrix parameters as well. 
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Fig. 3.5 Triangular Panel Data File Structure 

With these new parameters For the perspective transformation 
equations, the georectangular coordinates of the terrain 
grid are once again run through the perspective transfor- 
mation, and the new IA and JA coordinates are calculated for 
each point. These new IA and JA coordinates represent the 
transformed elevation paints as seen from the new observer 
location. The next step is to take the same three elevation 
paints that formed a triangular panel in the original 
transformation, map them into the synthesized image file 
using the new IA and JA coordinates, and then fill them in 
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with the assigned grey scale value. The mapping and filling 
process is explained in the next section. 

2 . Polygon Fill and Hidden Surface Elimination 

In the perspective transformation of the terrain 
grid into the synthesized image plane, some of the tri- 
angular panels became partially or entirely hidden by ether 
panels. To determine which pixels will be visible and 
therefore written to the synthesized image and which pixels 
are hidden, the z-buffer algorithm was used. 

When the new observer location is input into the 
program, the XL, YL , and 2L georectangular coordinates are 
tabulated. Using the georectangular coordinates previously 
calculated for each of the terrain grid elevation points, 
the distance or depth from the observer location to the 
elevation data paints can be calculated by the following 
equation 

Depth = Square root C C XL-X ) 2 + CYL-YD 2 -KZL-Z) 2 ) C3.143 

The depth of the pixels within a triangular panel 
were calculated by using the normalized plane equation that 
defines the plane of the triangular panel. The normalized 
plane equation in 3D space is given by 
ax + by + cz - -1 C3.155 

To solve for the coefficients a, b, and c, the three eleva- 
tion points that specify a triangular panel are used, and 
the x, y, and z are replaced with IA, JA, and the Depth of 
each elevation point. If the three paints are CIA1, JA1, 
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Depthl), (102, J02, Depth2), and (103, J03, Depth3), then in 
matrix form we have the fallowing 
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Salving far the coefficients a, b, and c we have CRef. 2, p. 
2001 
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The inverse af the elevation paint matrix is 
determined by calculating the adjoint, which is the trans- 
pose of the cofactor matrix, and multiplying each term by 
the reciprocal of the determinant CRef. 5, p. 0-151. The 
depth of each pixel within a tranformed triangular panel can 
new be determined by using the 10 and J0 of the pixel in the 
following equation. 

Depth - - Cl + a (10) + b (J0))/c (3.18) 

Each triangular panel will have different a, b, and c 
coefficients, therefore the pixels within each triangular 
panel will have a different depth than those from another 
panel . 

The z-buffer is a two dimensional array that is the 



same size as the synthesized image of 512 by 512 pixels. 



When the Iff and JA coordinates of a pixel is determined, the 
depth is calculated and then compared to any previously 
written depth in the z-buffer at that same coordinate 
location. If the depth is smaller, it means that the pixel 
is closer to the observer and would caver the previously 
written pixel. The depth of the new pixel will then replace 
that of the previous pixel. The assigned grey scale value, 
which is determined from the projected triangular panel in 
which the new pixel is contained, is written to the syn- 
thesized image at the calculated IA and JA coordinate 
location. If the new pixel depth is larger then the 
previous pixel depth, that means the new pixel is farther 
away from the observer and would be covered by the previous 
pixel. In this case no action is taken and the next pixel 
is processed. This procedure is continued until all the 
pixels within each translated triangular panel has been 
processed. CRef. 2, pp . 265-267] 

The IA and JA coordinates of the pixels in the 
synthesized image are calculated from each transformed 
triangular panel using an active edge list (AEL), J__Bucket, 
and frame buffer. The IA and JA screen coordinates of the 
three elevation points that define a triangular panel are 
used to generate the parameters of the AEL . Three lines are 
formed that connect each of the three translated elevation 
points and are identified as line 1, 2, and 3. A line is 
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defined between any two paints from which you can determine 
the I A associated with the maximum JA of the two points by 
IA__INCPT * C Ifi value of the maximum JA ) (3.19) 

How much I A changes for a one step change in JA is given by 

(3.20) 

which is the inverse slope of the line between the two 
points. The span of JA between two points is given by 

DELTA_JA - JA ( Point? ) - JACPoint 1) (3.21) 

For each line, the IA coordinate that corresponds to 
the maximum JA value of the line, the amount IA changes for 
each one unit step of JA, and the total span of JA is 
determined and stared into the AEL . After the line para- 
meters are placed into the AEL, the line identification 
number is put into the J Bucket , which is a one dimensional 
array of size 512, at the same location as the maximum JA of 
the line. In this way the J_Bucket acts as a pointer to the 
lines stored in the AEL. If a previous line number has 
already been written to that location, then the line ident- 
ification number of the new line is placed into the AEL and 
is linked to that line already residing in the J_,Bucket . An 
example is shown in Figure 3.6 where line #1 is referenced 
by the JBucket and line #2 is linked to line #1. The IA 
and JA coordinates of the pixel locations to be mapped to 
the synthesized image are contained within the three lines 



DELTA IA = 



IA(Point 2) - I A( Point 1) 



JA(Point 2) - JA(Point 1) 
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whose parameters are contained in the AEL . IF a line is 
horizontal, the IA and JA coordinates of that line are 
located between the other two lines and therefore does not 
need to be written to the AEL. CRef. 2, pp.75-7B3 

J_ Bucket AEL For Line 1 




AEL For Line 2 

Fig. 3.6 Active Edge List Storage 

The J_Bucket is scanned From its maximum to minimum 
coordinate value. IF the J_Bucket contains a line pointer, 
the parameters For that line are retrieved From the AEL as 
well as those of any line it is linked to. Since we have 
defined a triangle there will always be two lines to work 
with. From the parameters of the two lines, the IA coor- 
dinates that fall between them is calculated for each JA 
scan line. A scan line is all the columns CIA coordinates) 
along a particular row C JA coordinate) . As the JA coor- 
dinate is decremented by one, from maximum to minimum, the 
J__Bucket is checked to see if a new line has been added, and 
then the corresponding IA for each line is determined For 
that JA value. Those IA values and any that Fall between 
them are mapped to the Frame buffer. 
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The Frame buffer is a 512 by 512 array that is 
initialized to 0. As the IA and JA values are mapped into 
the frame buffer the coordinate locations matching the IA 
and JA values are changed From 0 to 1 . This process 
continues, and each time the JA scan line is decremented, 
the DELTA _J A parameters for each of the tuo lines are also 

decremented and the J Bucket is checked. If the J_Bucket 

contains another line pointer, then the line it references 
will replace the line whose DELTA _JA has decreased to 0. 

When finished, the frame buffer will have recorded the IA 
and JA coordinates for every pixel location within the 
translated triangular panel . The frame buffer is then 
scanned, and if a particular IA and JA coordinate location 
contains a value of 1, then the depth for a pixel located at 
those same coordinates is calculated and compared to 
previously tabulated depths in the z-buffer as explained 
earlier . 

This process is known as a polygon fill routine that 
utilizes the z-buffer algorithm for hidden surface elimina- 
tion. If any part of a triangular panel, after going 
through the perspective transformation, should map outside 
the synthesized image coordinate boundary, the entire panel 
is discarded and the next panel is processed . This is not a 
satisfactory solution and could have been corrected by 
implementing a clipping routine that would allow partial 
triangular panels to be mapped to the synthesized image. 
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However , due to time constraints, a clipping routine was not 
incorporated into the 30 transformation program. 

D. SUriMARY 

Using the artificial reference image from Figure 3.4 and 
an artificial terrain grid that mapped to the same image, 
two synthesized views were generated. The first synthesized 
view shown in Figure 3.7 was from an observation point 
located at 37° EE ’ 30” N. latitude, -1EE° 01’ 59” W. 
longitude, and an elevation of 110 meters. Figure 3.8 
exhibits the next synthesized view that was produced from an 
observation point of 37* E0 ’ 0” N. latitude while main- 
taining the same longitude and observer height as before. 
This simulates viewing the object from further away. As 
expected of a perspective view the object appears smaller. 
The near view also demonstrates the perspective relationship 
by the apparent taper from front to back. A third perspec- 
tive view is depicted in Figure 3.9 from close in and at a 
higher elevation. The observer location was from 37* S3’ 
30’’ N. Latitude, -1EE* 01’ 59’’ W . Longitude, and a height 
of 150 meters. This demonstrates the visual effect of not 
displaying partial triangular panels in the image and why a 
clipping routine is needed. 

The actual reference image is shown in Figure 3.10 from 
which the synthesized view displayed in Figure 3.11 was 
generated. In comparing the reference image to the syn- 
thesized view there is little resemblance. A closer 
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synthesized view is seen in Figure 3.12 which gives a 
partial representation of the reference image. The lack of 
detail in both shading and the depiction of features can be 
attributed to the poor resolution of the terrain data (ap- 
proximately 25 meter resolution) as compared to that of the 
reference image data Cl meter resolution). 

To improve the quality of the synthesized view the 
terrain data must be of a higher resolution. Having more 
elevation data points over a given area, the fewer number of 
reference image pixels one must collect and average for a 
triangular panel . The synthesized image will then maintain 
a closer approximation to the various shades of the refer- 
ence image. The physical shape of an object also suffers 
from poor terrain data resolution. The perspective trans- 
formation forms straight lines between translated elevation 
paints. This causes object distortion if the elevation 
points do not fall exactly along the boundaries of the 
object. As an example, if a square box in the reference 
image had only one elevation point defined on its surface at 
the center of the box, then the synthesized image can not 
reproduce the corners and edges of that box, and it would 
appear distorted. For these reasons the closer the resolu- 
tion of the terrain data to the resolution of the reference 
image, the closer the synthesized view resembles the refer- 
ence image in shading and shape. 
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Fig. 3.7 Artificial Synthesized Image 1 




Fig. 3.9 Artificial Synthesized Image 2 
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Fig. 3.9 Artificial Synthesized Image 3 




Fig. 3.10 Reference Image 
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Fig 3.11 Synthesized Image 1 



Fig. 3.12 Synthesized Image 2 
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The 3D transformation program was developed to allow 
tracing the Flow of data easily. Much of the data was 
written to files so that it could be printed out and 
studied. This method of data storage required certain files 
to be read many times as the synthesized image was gener- 
ated. The program execution would have been faster if the 
data had been stored in arrays that are easily passed 
between the various modules. Another consideration that 
would increase speed would be to decrease or eliminate the 
interactive input required. This could be accomplished by 
extracting the desired terrain data into a file and sepa- 
rating the associated reference image before program 
execution. This would require longer set up time but would 
decrease the amount of data manipulation required by the 
program and increase its overall speed. 
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IU . CONCLUSIONS 



A . GENERAL 

The 3D computer image transformation from a photographic 
image mas a difficult task to achieve satisfactorily. The 
main objective mas to develop a program that takes a grey 
scale photographic image, a set of elevation data points 
defined over that image, and generates a rotated synthesized 
perspective view. This goal was realized, but some areas 
still need improvement. The quality of the synthesized view 
is judged on how well the grey scale values of pixels match 
those of the original image and how closely the translated 
objects resemble the desired structure. 

The resolution of the terrain model as compared to the 
resolution of the reference image data, extensively affects 
the synthesized image quality. Depending on the contents of 
the reference image, the required resolution of the terrain 
model will vary. If the image is of open country, the 
distance between elevation points may be large and the 
result may not suffer unacceptable degradation in the syn- 
thesized view. If the image is of a city that has many 
small distinct objects such as buildings, then the resolu- 
tion of the elevation points must be higher. Given a set of 
elevation data points, the question to be resolved is how to 
make the synthesized pixels relate to the reference image 
better, and thus improve the quality of the synthesized 
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views? This question and alternative ways to improve speed 
and program flexibility will be discussed. 

B. GREY SCALE CORRELATION 

There are a few possible methods to improve the shading 
of the synthesized view to better match that of the original 
image. The translated triangular panels form very distinc- 
tive lines or boundaries between areas of different shades. 
It may be desirable to blend these boundaries by sampling 
the pixels along both sides, replacing those pixels with an 
averaged value. Another possibility would be to implement a 
Gouraud or Phong shading algoritm CRef. 3, pp . 323-330J . 

This would help smooth the shade transition across the 
boundary and result in a more smooth appearance. 

Another method to improve grey scale correlation would 
be to develop a way to divide the triangular panels into 
smaller triangular areas before referencing the original 
image. By having smaller triangular panels, smaller areas 
of the reference image are sampled resulting in a better 
approximation in the synthesized view. 

Only the left image of a stereo photographic pair of 
images was used in this study. The right image may contain 
grey scale information of surfaces hidden in the left image 
view cr vice versa. By sampling both left and right images 
and complementing them, some ambiguities may be resolved. 
These suggestions will increase the number of calculations 
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required to generate the synthesized view but may improve 
the quality of the synthesized view. 

C. PROGRAM SPEED AND FLEXIBILITY 

Although speed was not a prime consideration in this 
study, it would be desirable to generate synthesized views 
as quickly as possible. To determine which subroutines 
consumed the most CPU time the program was monitored while 
running the artificial reference data set. Table 1 shows 
the results of the CPU time used and the percentage of the 
total time for each subroutine called by the main program. 

If a subroutine calls another subroutine that time is 
included in the calling routines CPU time. The subroutines 
that required interactive input were not measured. There- 
fore, the total CPU time cannot be measured precisely. 
However, the results obtained represent reasonable estimates 
and demonstrate where the focus should be on improving 
overall speed. 

Some suggestions have already been discussed on how to 
improve the speed of the program, but other methods also 
exist. Preconditioning the input data to eleminate interac- 
tive input and using arrays instead of files were two 
methods already presented. The z-buffer is accessed several 
times during synthesized view generation. If the z-buffer 
could be implemented in hardware, the time required for 
processing of polygon fill and hidden surface elimination 
routines would be improved. 
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TABLE 1 



CPU TinE CONSUMPTION CIN SECONDS) 



SUBROUTINE 


CPU T 


TER _ CROP 


14 . 


READIMAGE 


3. 


TER INTRP 


0 . 


REAL EL 


0 . 


TER_ DMS 


0 . 


IM REF I J 


0 . 


I M__REFAUG 


1 . 


AFFIN 


0 . 


NEW_I J 


0 . 


NODE _DEPTH 


0 . 


FILL 


134 . 


TOTAL 


155 . 



I ME 


PERCENTAGE 


11 


3.051 


23 


2.072 


04 


0.026 


07 


0.045 


63 


0.443 


33 


0.250 


10 


0.706 


04 


0.026 


00 


0.513 


63 


0.404 


80 


86.467 



3 100.0 



The Frame buffer is used and accessed in two different 
areas of the program for each triangular panel translated. 

It may be possible to eliminate this buffer by processing 
each pixel as its IA and JA coordinates are determined, 
instead of storing that information into the frame buffer 
for later processing. Other polygon fill routines such as 
seed fill algorithms or using fence registers may be faster 
than the edge fill routine used in this program CRef. 2, pp . 
30-BBJ . 

For program flexibility it uould be desirable to be able 
to adjust the image plane of the synthesized vieu to any 
desired angle. The program limits the image plane to a 
northern direction. To make the synthesized image plane 
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adjustable mould require developing a method for rotating 
the image plane coordinates in terms of the georectangular 
coordinates. This mould improve the 3D transformation 
program and extend its usefulness. Another program improve- 
ment mould be the incorporation of a clipping routine to 
improve the appearance of partial synthesized viems as 
discussed in previous sections. 

The ideas used to develop this program mere contrived 
from fundamental concepts . Plany areas mere discussed that 
could improve or enhance the basic implementation of the 
program. The results that mere obtained are encouraging and 
could easily be used as the basis for further study. 
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*> o — CO 



PPPENEIX A 

PRCcPfi^ S'jnnARV 

PFCCP^h TP^ , \ I 3 3 
3. * crs ^s^Fdirrrscl 

Takes an i mage and slsvaticn data fils ard extracts 

ratsfcrmatitn . Accepts interactive 



— sa ^ 

' — *• — y - w - w - 



iput cf a *"eui observer location and generates a synthesized 



ISO. w ^ “ioi.* 






to various subrouting meddles 



te mai^ part cf the program 



e . 



Passes parameters betuee° subroutines , 

iriPGES file that contains the synthesized vie., 
Calling routines 
hone . 

Called routines 

I hPL’T : Obtains the interactive input cf 

the elevation a^d image file names 
and desired reus and ool_n n s ^sed 
to s:. tract slevatio^ data ooi^ts 
5 non t w e reference elevatio" gr^d 
E -tracts the desired elevation 
coirts cf the area to be tra^s- 



T E P CFCP : 






READ I HAGE : 


Reads the reference image intc an 
array called IhACE. 


yrp t p.j-r oo . 


This rc^t^re cciiacts and averages 
the extracted elevation pants and 
assigns that va Ice tc any unkmciun 
eievaticn data pcnts, 


REAL EL: 


L rites the extracted elevation 
points into a file called ZFIL.DAT 
as real values m meters abcve sea 
level . 


TEP DCS: 


Determines the latitude and lcn- 
gitude cf the eievaticn points then 
converts them tc gecrectangular 
coordinates and stores them m a 
file called XYZ . DAT . 


T M per T ' r 


Converts the georectangular coor- 
dinates to the reference image I.A 
and JA screen coordinates. 


ir REFAUG: 


Constructs the file MODE. CAT that 
contains the three elevation points 
and grey scale value that make uc a 
triangular panel . 


AFFIN: 


Determines the affine transform 
coefficients utilized i n the 
transform cf the synthesized image 



CQ 



npr r np 



rjE’jJ IJ: 



kj n ^ r~ n o t u . 



plane coordinates to screen 
coordinates . 

Accepts the interactive input cf 



the latit 


, — J -2 


longitude, and hei 


in meters 


abo' 


e sea level of the 



desired observer location. 

Computes the reu IP and JP screen 
coordinates of the elevation points 
for the synthesized view. 

Calculates the distance from each 
elevation point to the observer 
looat ion . 



FILL: Determines the hidden surfaces and 

generates the synthesized image in 
the file I r-PGES . DPT . 

Routine parameters 

I EM DM : The number cf rows cf the extracted 

elevation data points. 

I El T i 0 hi : The number cf columns of the 

extracted elevation data points. 



E . SUBROUTINE INPUT 

a. Functions performed. 

Pccepts the interactive input of the elevaticr a~d 
^ mage data files, as cell as the information needed to 
the desired ele v a 1 1 c n points. 



e < crac t 



Input (interactive) 



IRCUJ: 


Nmimu n rou of the e 
desired . 


levation area 


LRCliJ : 


Maximum rou of the s 
desired . 


levaticn area 


I CCL : 


Minimum column value 
elevation area . 


cf the desired 


LCCL : 


Maximum column value 
elevation area . 


of the desired 


ELFILE : 


The reference slevat 


ion file name . 


I riF I LE : 


The reference image 


file name . 


I FRAPE : 


The I Frame value of 
image . 


the reference 


JF P A?1E . 

GutpLjt 


The J Frame value of 
image . 


the reference 


Same as the 


interactive input . 





Calling routines 
TRGN 3 3. 

Called routines 
N one . 

Routine parameters 
None . 



3. EL! SR CUT I NE TER CRGP 



a . 



r u^ct ions 



7^C= P S 



H 
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Reads the original terrain grid and constructs a 
smaller grid of the desired elevation points m the array 



IELEU2 . 

b . I nput 
T PC'J • 

LRGU; : 



LCGL : 



ELFIL: 



Minimum rcu of the desired 
elevation area, 
flaximum rou of the desired 
elevation area. 

n i n i m o m column value of the desired 
elevation area . 

Maximum column value of the desired 
elevation area . 

The file containing the reference 
elevation data. 



o . Output 

lELE'JE: An array containing the extracted 

elevation points. 

d. Galling routines 
7PAN 3 D. 



Galled routines. 



More . 



Routine parameters 

IELEU1: An array containing the entire 

elevation data set. 



and 



Counters . 



p 



4. SUBROUTINE HEAD I MAGE 

a. Functions performed 

Reads the reference image into an array called 

image . 



TM|T T T . 



output 



IMAGE: 



The file containing the reference 
image pixel grey scale data. 

An array containing the pixel grey 
scale values of the reference 
image . 



Galling routines 
TRAM 3 D . 

Galled routines 
None . 

Routine parameter 
IE and I G : 



counts: 



5. SUBROUTINE TER INTER 

a. F -notions performed 

Determines the average elevation over the extracted 
elevation area and assigns that elevation to any unknown 
elevation points, 
b . Input 

I ENCM : "lumber of extracted elevation reus 

IENGM: Number of extracted elevation 

co limns . 



cm 



IELEU2: 



°r array cf the extracted elevat 
data . 



wotpUt 

I ELEUE : 



^rr3 a containing 


t ^ e 


e ' t r 


acted 


elevation data 


sny 


un knoun data 


points have bse n 


ass 


•> m a — 1 1 

- 3 ' — — 


the 


average elevat i a 


n of 


the 


area . 



Ea ill ng 



rcjtmes 



TEAM 3 D . 

2a 1 led rr jt i^ss 



here : 

Routine parameters 
N and r 1 : Counters . 

NEXT: Count cf the number cf elevation 

points processed. 

I EL : Summation cf the elevations . 

[ ~UG : The avsrags elevation cf the area. 



E. 5LERCLT I NES REAL EL 

a. Functions performed. 

Creates a real file called ZEIL.GAT cf the extract ad 
c 1 2 a*~'0 T ~ points, 
b . I rput 

I EMDf'j : ^he number of extracted elevatio^ 

reus . 

I ELCL : The number cf e.; treated elevation 



oo * jr 



E4 



IELEU2: An array cf the extracted elevation 

data pcirts . 

Output 

ZFIL.DAT. A file rcntai-irg the real value cf 

the extracted elevation data 
p C 1 t s . 

falling ri_ti n ss 

T?PN 3 □ , 

failed routines 
None . 

Routine parameters 

AELEU : — real array of the extracted 

elevation points. 



SUSRCTJT I r JE TER DrS 



a. Functions performed 

fen i' arts each extracted elevation data pci^t to its 
frS equivalent. It uses the fact that each elevation point 
represents a cne second change in latitude cr longitude f r cm 
the next point, starting f r cm the Southwest corner reference 



ele 5 


■aticn located 


at 37 


° EE 


' 4” " M. Latitude and -122* C5 


w 1 w 


W. Longitude. 


The 


□ris 


is converted to X, Y and 2 



gscrectangular ccordinates and stcred i n the XYZ.fAT file, 
h. Input 



The minimum row cf the extracted 
elevation data prints. 
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LEQU: 



The .Tax i mum rou cf extracted eleva- 
tion data points. 



-J 

□ 

u 

t — 1 


The 


minimum column 


value 


cf the 




extr 


acted 


elevation 


data 


points 




The 


maximum column 


value 


c f the 




extr 


acted 


elevation 


data 


points 


IE-iDr : 


The 


total 


number cf 


reus 


cf the 




extr 


acted 


elevation 


data 


points 


ZFIL.DAT: 


The 


file c 


cntainirg 


the r 


sal v a 



of the extracted elevation data 
points . 

Output 

XY2 . CPT : P file containing the gecrec- 

tangular coordinates cf the ex- 
tracted elevation data points. 

Tailing routines 

TRPN 3 D. 

Tailed routines 
CHS2XYZ . 

Routine parameters 

IL and JL: Counters. 

X, V, and Z: Geor ectangui ar coordinates of an 

elevation point . 

L^TT . LPT M , and 

T-TG: The latitude i n degrees, minutes, 

a^d seconds of an elevation point. 



EG 



LGND, LONH, 

and LONS: The longitude in degrees, minutes, 

and seconds cF an elevation point. 



HEIGH 


-r 




The 


real elevaticn in meters 


above 








sea 


level cf an elevaticn pci 


* 4- 


PLAT” 


and 


PLATS: 


The 


reference latitude in mm 


Li t BS 








and 


seccnds cf the reference 


eleva- 








tier 

i 


pemt lacated at rcu 1, 


c c i j m n 


PLONN 


AND 


PLONS : 


Tlnn 
x J : tr 


reference longitude in mi 


nutes 








and 


seccnds cf the reference 


eleva- 








t izn 


pcint lecatsd at rcu 1, 


cclumn 



B. SUBROUTINES in. REF I J 

a. Functions performed 

Calculates the I and J reference image coordinates 
and converts them to the IA and JA screen coordinates for 
each extracted elevation data point and puts them into 
arra 3 s IA, and w T A. 
t . Input 



I FRAME: 



JFRAME: 



The I Frame value of t K e reference 
image . 

' T 'Ha J" frame VSiuS IjP t h S csfsrs PC3 
ir-ags . 

The number cf reus in the extracted 



e i s a 



ion 



data . 



IENDN: 



The number o f columns in the ex- 
tracted elevation data. 

XYZ . OPT : Tata file containing the g e c r e c - 

tangjiar coordinates of the 
extracted elevation data points. 

Cut put 

IA: An array containing the IP screen 

coordinate values of the extracted 
elevation data points. 

JA: An array containing the JA screen 

coordinate values of the extracted 
elevation data points. 



Calling routines 








TFAN 3 □ . 








Called routines 








PRCJECT XYSIJ. 








Routine parameter 


•s 






GHEE A : 


Rotation cf 


the 


original image 




plane about 


:>•: — a ; 


x i s i n radians. 


PHI : 


Rctaticn cf 


the 


original image 




plane about 


the 


g - a x i s m radians 


KAPPA: 


Rctaticn cf 


the 


original image 




pla^e about 


the 


z-axis m radians 


XC and YD: 


Cffset cf t h 


e principal pcint frc: 




the IPP. 







E0 



and Z I : 



The 


gsorectangular 


CO 


rt 

tr 

0) 


observer locat 


ion 



ai , 


A2, 


P 1 

w — , 


52 , 




Cl , 


and 


CZ : 




The given affine transform para 



meters for the transform of the 
original image plane coordinates to 
I and J original image screen 
coordinates . 



FGCL'S : 


The camera focal 


length 


• 


: 


The -navis value 


cf the 


elevation 




point m image p 


lane cc 


ordinates . 


Y I : 


The Y-axis value 


r-> P h e 


elevation 



point m image plane coordinates. 
I and J: The original reference image 

screen coordinates. 



a t 5 U E ROU T I N E nn S 2 X V 2 

a. Functions performed 

Converts the DHS latitude a^d longitude data to X, 
Y, ^nd Z gecrectangular coordinates. 



T p i +• 

- 1 j- — - 



LAIC, LPin, and 
L^TS : 



The latitude in degrees, ni^utss 
and seconds cf the e x t r a c t s d 
elevation data points. 



( r I M 



cq 



LCMS : The longitude in degress, m i n u t e s 

and seoo^ds of the extracted 
elevation data points . 

HEIGHT: The height i n meters above sea 

level cf the extracted slevatio~ 
data points . 



X, V, and Z: The gecrectangular coordinates of 

t, h e e . t r a o t e d e 1 s % ■ a 1 1 o n c * ^ s . 

Tailing routines 

~r p* P rp p-1 CQ HOC T r> £ 

Called routines 
Cicre . 

Fouitine parameters 

PHI: Angle m radians frcm tte 2-a:;:s rf 

the gsorsrtanguiar rerrdi rate 



L~r ; QA : 



N , E . SQUARE 
and A : 



P^n ▼ n i 

P I : 



system . 

Angle m radians from the X-axis of 
the georeot angular coordinates 
system . 

Given parameters used m calcu- 
lating the georecta ngular coor- 
dinates of the elevation points. 
Number of radians per degree. 

Number of radians i ~ a half enrols. 



lumber 



af degrees in a half 



Cl : 



e . 



CE: Number cf minutes in a degree. 

C3-. Number af seccrds in a degree. 

1C . SL'EPGUT I HE P P □ J E E T 

a. Functions performed 

Converts the X, Y, sot 2 gecrectangular coordinates 
of the extracted elevation data points to the x- image ant 
g-ioags coordinates of the mage plane, 
b . I npu t 

Of'lESfP: x*a>;:5 rotation cf the image plane 

in radians . 

PHI: y-axis rotation cf the image plane 

in radians. 

KPPPS: z-axis rotation cf the image plane 

in radians . 

XG and YQ : The offset cf the principal point 

f ram the I PP . 

XI, Yi, and 21: The gecrectangular cccrdmates cf 

the observer location. 

X, Y, and 2: The gecrectangular coordinates cf 

the extracted elevation data 
pc m ts . 

FGCUS: The camera focal length. 






XIHA: 



The image plane x-axis coordinate 



value of the extracted elevation 
data points. 

Y I H A : The image plans y-axis cccrdinate 

value st the extracted elevation 
paints . 

Tailing routines 

in pefij. 

Tailed routines 

None . 



Routine parameters 



fill , 




n 13 


M21 , 


r-jgp 


riE3 


n3 1 , 


m3E 


M 33 


r-jc-Mnr 


i . 





The parameters of the n matrix. 

The denominator of the transforma- 
tion equations. 



RGL'TI.NE XYEIJ 
Functions performed 

Ton verts the image plane coordinates of the ex- 
elevation data points to the I and J original image 
ates . 

I nput 



XI HA : 



The image plane x-axis coordinate 
values of the extracted elevation 



data points . 



Y I n A : 



The image plane y-axis coordinate 



values of the extracted elevation 
data pc in ts . 

Al, AS, SI, ES, 

Cl, CB: The affine transform parameters to 

map the image plane coordinates to 
the reference image I and J 
coordinates . 

c . Output 

I and J : The reference image coordinate 

values for the extracted elevation 
data points . 

d. Calling routines 

in REF I J , NEW IJ. 

s . Called routines 

None . 

f. Routine parameters 

TENON: The denominator cf the affire 

transf crm . 



12. SUBROUTINES I N-REFAUS 
a. Functions performed 

Constructs the file NOTE. CAT that contains the 
elevation points and reference grey scale value that make up 
a triangular panel. 



Input 






I A and 



im^se : 

I ENOn : 

I ENON : 

Output 
NODE . EA 



Oal : mg 
TRAN 3 
Called 
None . 
Routine 
SLOPE : 



JA: The arrays that contain the re- 

ference image screen coordinates of 
the extracted elevation data 
points . 

The array that ccntains the 
reference image grey scale values. 
The number of rows m the extracted 
elevation data. 

The number of columns m the 
extracted elevation data. 

T: The file that contains the eleva- 

tion points that make up a tri- 
angular panel and its associated 
reference grey scale value. 

routines 
0 . 

routines 

parameters 

The slope cf the diagonal line that 
separates the rectangle, defined by 
four elevation data points, into 
two triangular panels. 

The Y intercept of the diagonal 
line. 
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IV, n, and N: 
I GREY 1 : 



I3PEYS: 



L, LI, and IP: 
MGOE A, NODE 5, 
NODE C, and 
NCEE D: 



I COUNT 1 : 



ITOT1 : 



I COUNTS : 



Counters . 

The averaged reference grey scale 
value of the triangular panel belcw 
the diagonal line. 

The averaged referenced grey scale 
value cf the triangular panel above 
the diagonal line. 

Counters . 



The numerical designation cf the 
four elevation points that make up 
the rectangle that is divided into 
two triangular panels. 

The count of the number of pixels 
averaged in the triangular panel 
below the diagonal line. 

The summation of the pixel grey 
scale values averaged m the 
triangular panel below the diagonal 
line. 

The count of the number of pixels 
averaged in the triangular panel 
above the diagonal line. 

The summation of the pixel grey 
scale values averaged m the 



-7C; 



triangular panel above the diagonal 
line. 



13. SUBROUTINE EXT ORIEN 

a. Functions performed 

Determines the fl matrix parameters used m the 
rotation of the image plane for the synthesized view. This 
program defines the three image plane coordinates m terms 
of georeot angu 1 ar coordinates. 

b. Input 
IENDM: 

I ENDN 

o . Output 

r 1 1 , ms, n 1 3 : 

nai, ri22, ri33: 
r3i, ri3B , ri33 : 

d. Calling routines 
New I J . 

e. Called routines 
None . 

f. Routine parameters 



The number of rows in the extracted 
elevation data points. 

The number of columns in the 
extracted elevation data points. 

First row coefficients of the 
transformation matrix. 

Second row coefficients of the 
transformation matrix . 

Third row coefficients of the 
transformation matrix. 
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HAGN X: 


The magnitude c F the synthesized 
vieu image plane x-axis. 


flAGN Y: 


The magnitude cF the synthesized 
vi su image plane y-axis. 


HAGM 2: 


The magnitude cf the synthesized 
view image plane z-axis. 


X, Y, and 2: 


Georectangular ccordinates cF the 
extracted elevation data points. 


X CORD: 


Array that stores the georec- 
tangular X coordinates cF the 
extracted elevation data points. 


Y CCRD: 


Array that stares the georsc- 
tangular V coordinates oF the 
extracted elevation points. 


2 CCRD: 


Array that stcres the georec- 
tangular Z ccordinates cF the 
extracted elevation data points. 



UECX, X UECY 



and X UECZ: 


The X, V, and Z georectangular 
cocrdi nates oF the synthesized 




image plane x-axis. 



UECX , Y UECY 



and Y UECZ: 


The X, Y, and Z gecrectangoiar 
coordinates of the synthesized 




image plane y-axis. 
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2 UECX, 2 UECY , 

and 2 UEC2 : The X, Y, and 2 georec tangul ar 

coordinates of the synthesized 
image plane z-axis. 

ITCT: The total number of extracted 

elevation paints. 

IP and IP: Counters. 

14. SUBROUTINE AFFIN 

a. Functions perfcrmed 

Assigns or calculates the coefficients to be 
utilized m the affine transform from image coordinates to 
screen coordinates of the synthesized view. 



b . 


I npjt 








None . 

n, .i-n, i- 








-UUPL U 

PI, AP, 


51 , 52 






and C 1 , 


C2: 


The affine transform coefficients 



for the synthesized view image 
plane coordinates to screen 
coordinates . 



d. Calling routines 
TRAN 3 D . 

e. Called routines 
None 

f. Routine parameters 



p 



X 1 11 A 


mx : 


Assigned maximum 
coordinate . 


image plane x 


VI HA 


hi AX : 


Assigned maximum 
coordinate . 


image plane y 


T M/S V 




Assigned maximum 
coordinate . 


IA screen 


j mx 




Assigned maximum 
coordinate . 


JA screen 



15. SGERDUT I PJE DBS LOC 

a. Functions performed 

Calculates the neu observer location georectanguiar 
coordinates from the interactive input of the desired 
latitude, a^d longitude, and height. This routine also 
assigns the focal length for the synthesized view, 
b . Input ('interactive) 

LATH and LATS: The minutes and seconds of the 

latitude of the neu observer 
location . 

LCN M and LONS: The minutes and seconds cf the 

longitude of the neu observer 
location . 

HEIGHT: T H e altitude in meters above sea 



level f 



the neu observer lcca- 



XI, VI, and 21: The gecrectangu 1 ar coordinates of 

new observer location. 

FOCUS: Assigned focal length for the 

synthesized view . 

d. Tailing routines 
TRAN 3 D . 

e. Tailed routines 
DMS2XV2 . 



Routine parameters 



None 



IB. SUBROUTINE NEW IJ 



a. Functions performed 

Calculates the new IA, and JA screen coordinates of 
the extracted elevation data points using the new observer 
location data . 

b. Input 



XI, VI and 21: Gecrectangular coordinates of the 

new observer location. 

FOCUS: Assigned focal length for the 

synthesized view. 



Al, AS, Bl, BS , 
and Cl, CE: 

I ENON : 



The affine transform coefficients. 
The number of reus of the extracted 
elevation data points. 



SO 



IENDN: 



The number of columns of the 



extracted elevation data pcmts. 



Output 



IA: 



The array containing the elevation 



data IA screen cccrdinates cf the 



synthesized view . 



JA : 



The array containing the elevation 



data JA screen cocrd mates cf the 



synthesized view . 



d, Tailing routines 
TRAN 3 D. 

e. Called routines 
XY2IJ, EXT QRIEN. 

F. Routine parameters 
HU, HIE, and 

n 1 3 : First rout coefficients cf the 



transformation matrix. 




H33: 



Second rou coefficients of the 



transformation mat r l 



M31, M32, and 



Third rou coefficients cf the 



transformation matrix. 




The offset of the principal point 



to the I PP . 



Image plane x coordinate. 
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YIHA: 



Image plane y coordinate, 



DENOn : 



Denominator of the transformation 
matrix . 



ITGT: Total number of extracted elevation 

data paints . 

I R : Gaunter . 

1“\ SUBROUTINE NCDE EPTH 

a. Functions performed 

Calculates the distance from each transformed 
elevation data point to the neui observer location. 

b. Input 



XI, Y1 , and 21: The gear ec tangu 1 ar coordinates of 

the neuj observer location. 

I ENOn : The number of rows m the extracted 

elevation data points. 

IENDN: The number of columns m the 

extracted elevation data points. 

Output 

DEPTH: Array of the distances from the 

transformed elevation data points 
to the new observer location. 

Calling routines 



TRAN 3 D. 

e. Called routines 



None . 

Routine parameters 
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IP: 


Haunter . 


I TGT : 


Total number of extracted elevation 
data points. 


»q C 'pD n ^ rM !T ^ T * t 





a. Functions performed 



Determines the 


hidden surfaces using a z-buffer 


algorithm. The depth : 


is compared for each pixel at a 


specified screen coord: 


in ate to determine if it is written to 



the synthesized image 



b. Input 




I P : 


Prray containing the IP screen 
coordinates of the transformed 
elevation data points. 


jp : 


Prray containing the JP screen 
coordinates of the transformed 
elevation data points. 


CTPTH : 


Prray of the distance from the 
transformed elevation data points 
to the new observer location. 


I HPGE : 


Prray used to construct the 
synthesized image . 


lENDr : 


The number of rows of the extracted 
elevation data points. 


I ENGN : 


The number of columns of the 
extracted elevation data points. 



S3 



XYZ.DAT: 



File Qf the gearectangular cacr- 



dinates nf the elevation data 
points . 

Output 

IPPGES : P file containing the synthesized 

image . 

Calling routines 
TRPN 3 D . 

Called routines 
FRPHE F I L . 



Routine parameters 
Z EUFF: The z-huffer 



y i 



i 1 2i X2 
YE Z P V 3 Y3 

I w , , - — - 1 * 1 



and Z3: 



1PTU - 



The screen coordinates and depth of 
the three elevation points that 
define a triangular panel. 

The distance of a pixel uithm a 
triangular panel to the new 
observer location . 



r 7 i 


C 12 , 


Cl 3 , 




C21 , 


C22 , 


C23 , 




C31, 


C32, 


C33 : 


The cofactors of the plane equation 



matrix. 



GET: 



The determinate of the plane 
equation matrix. 
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A COEF, S COEF, 



and C.COEF: 


The coefficients of the plane 
equation . 


I GREY: 


Grey scale value of a triangular 
panel . 


NODE A, NCDE 


p 


NODE C : 


Numerical designation of the three 
elevation data points that make a 
triangular panel . 


T M T V] . 


ninimum IA screen coordinate of 
the three elevation points. 


I PI AX : 


Maximum IA screen coordinate of the 
three elevation paints. 


J MIN: 


Minimum JA screen coordinate of the 
three elevation points. 


J MAX : 


Maximum JA screen coordinate of 
the three elevation points. 


FRAME : 


The frame buffer. 



IP, I , J 



K, and L: 


Counters . 


I PLANES : 


Total number of triangular panels 
constructed from the extracted 
elevation data points. 



PC 



19. SUBROUTINE FRAME FIL 



a. Functions performed 

Constructs an edge list from the three transfcrmed 
elevation points. These edges form the boundaries for a 
polygon fill routine. The pixels that are determined to 
fall within the transformed triangle are marked m the frame 
buffer . 



H 



d . 



I npu t 



NODE A, NODE B, 

and NODE 0: The numerical designation of the 

three elevation data points that 
make a triangular panel . 

I A *. Array that contains the IA screen 

coordinates of the transfcrmed 
elevation data points. 

JA: Array that contains the JA screen 

coordinates of the transformed 
elevation data points. 



J MAX: Maximum JA screen coordinate of the 

three elevation data points. 

J MIN: Minimum JA screen coordinate of the 

three elevation data points. 



Output 

FRAME: The frame buffer array. 

Calling routine 
FILL . 



8E 



Called routine 



e . 



None . 

Routine parameters 



I INCPT: 


The matching IP coordinate cf the 
maximum JP point of a line. 


DELTA I : 


The amount IP changes for a one 
step change m JP . 


DELTA J: 


The total JP span cf a line. 


AEL: 


Prray containing the parameters of 
the three lines cf transformed 




triangular panel. 



XNDDE1 and 



XNODE 2: 


Designates the number of IA coor 


dmatss between two 


lines For a given 



J P 



DX: 


Indicator used to determine the 
direction in which the IP coor- 
dinates are counted. 


NODE: 


Prray containing the numerical 
designation of the three elevation 
data paints . 



NDDE1, and 



N0DE2: 


Designators for determining the 
elevation point that contains the 
highest JP value between two 




points . 
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N HIGH and 



N ..LOW: 



IT, I U , IP 
and IS: 

J 5UGKET: 



II and IS: 



ICNT and LCNT: 



Used with MOOE1 and NODES for 
determining the highest JP value 
between two points. 

Counters . 

^rray that contains the line number 
designators referencing the PEL. 
Used to determine the number of IP 
coordinates to be written to the 
frame buff er . 

Counters for the J EUCKET . 



SB 
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APPENDIX B 



3D -TRANS FORMAT I ON PROGRAM LISTING 
PROGRAM TRAN_3_D 

THIS PROGRAM TAKES AN IMAGE AND ELEVATION FILE AND 
CONSTRUCTSTHE REFERENCE IMAGE AND ELEVATION FILE. 
FROM THESE FILES A SYNTHESIZED IMAGE IS PRODUCED. 

CHARACTER ELFILE*13, IMFILE*13 
BYTE IMAGE( 512,512) 

INTEGER IELEV2( 50,50) , IA( 2500) ,JA( 2500) 

I NTEGER I ROW , LROW , I COL , LCOL , I ENDN , I ENDM , 

INTEGER I FRAME , JFRAME 

REAL XI, Y1,Z1, FOCUS ,DEPTH( 2500) 

REAL A1 , A2 , B1 , B2 , Cl , C2 , OMEGA, PHI , KAPPA 

CALL INPUT( I ROW, LROW, ICOL, LCOL, ELFILE, IMFILE, 

I FRAME, JFRAME) 

OPEN( UNIT=1 , FI LE=ELFILE , STATUS= ' OLD 1 ) 

OPEN( UNIT=2, FI LE=' ZFIL.DAT' , STATUS= ' NEW' , 

ACCESS=' DIRECT' , RECORDSIZE=128 , MAXREC=512 ) 
OPEN(UNIT=3,FILE='XYZ. DAT' , STATUS= ' NEW ' , 

ACCESS= ' SEQUENTIAL ' , FORM= ' FORMATTED ' ) 

OPEN( UNIT=4 , FILE=IMFILE , STATUS= ' OLD ' , ACCESS= ' DIRECT ' , 
RECORDSIZE=128,MAXREC=512 ) 

OPEN( UNIT=20 , FI LE= ' NODE. DAT' , STATUS= ' NEW ' , 

ACCESS=' DIRECT' , RECORDS I ZE=128, MAXREC=512 ) 

OPEN( UNIT=21 , FILE=' IMAGES. DAT' , STATUS='NEW' , 

ACCESS= ' DIRECT ' , RECORDSIZE=128 , MAXREC=512 ) 

CALL TER_CROP( IROW, LROW, ICOL, LCOL, IELEV2 ) 

I ENDN=LROW- I ROW+ 1 
I ENDM=LCOL- ICOL+ 1 
CALL READ I MAGE ( IMAGE) 

CALL TER_INTRP( I ENDN, I ENDM, IELEV2 ) 

CALL REAL_EL( I ENDN, I ENDM, IELEV2 ) 

CALL TER_DMS( IROW, LROW, ICOL, LCOL, I ENDM) 

CALL IM_REFI J( I A, JA, I FRAME , JFRAME , I ENDM, I ENDN) 

CALL IM_REFAVG( I A, JA, IMAGE, I ENDM, I ENDN) 

CALL AFFIN( Al, A2,B1,B2,C1,C2) 

CALL 0BS_L0C(X1,Y1,Z1, FOCUS) 

CALL NEW_ I J ( XI , Y1 , Z1 , FOCUS , Al , A2 , B1 , B2 , Cl , C2 , 

I ENDM, I ENDN, IA, JA) 

CALL NODE_DPTH( XI, Y1,Z1, DEPTH, I ENDM, I ENDN) 

CALL FILL( IA, JA, DEPTH, IMAGE, I ENDM, I ENDN) 

CLOSE( 1) 

CL0SE( 2 ) 

CLOSE( 3 ) 

CLOSE( 4) 

CLOSE( 20) 
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c 



c 

c 

c 

c 

c 

c 

c 

c 

c 



c 



35 



45 



55 



C 

C 

c 

c 

c 

c 

c 



c 



CLOSE( 21) 

END 

********************************************************** 
SUBROUTINE INPUT( IROW, LROW, ICOL, LCOL, ELFILE, IMFILE, 

I FRAME , JFRAME ) 



IROW : INITIAL ROW OF DESIRED AREA 

LROW : LAST ROW OF DESIRED AREA 

ICOL : INITIAL COLUMN OF DESIRED AREA 

LCOL : LAST COLUMN OF DESIRED AREA 

ELFILE : ELEVATION DATA FILE NAME 
IMFILE : IMAGE DATA FILE NAME 

I FRAME : I FRAME NUMBER 

JFRAME : J FRAME NUMBER 

CHARACTER ELFILE*13, IMFILE*13 
INTEGER IROW, LROW, ICOL, LCOL, I FRAME , JFRAME 



WRITE( 6, *) ' INPUT 
WR I TE ( 6, * ) 'ENTER 
READ( 5, 35) IROW 
WRITE( 6, *) 'ENTER 
READ( 5,35) LROW 
WRITE( 6, *) 'ENTER 
READ( 5, 35) ICOL 
WRITE( 6,*) 'ENTER 
READ( 5,35) LCOL 
FORMAT ( 13) 

WRITE( 6, * ) ' INPUT 
READ( 5,45) ELFILE 
WRITE( 6, *) ' INPUT 
READ( 5,45) IMFILE 
FORMAT( A13 ) 
WRITE( 6, *) ' INPUT 
READ( 5,55) I FRAME 
WRI TE( 6, * ) ' INPUT 
READ( 5,55) JFRAME 
FORMAT( 14) 

WRITE( 6, * ) ' ***** 

RETURN 

END 



ELEVATION AREA DESIRED' 

MINIMUM ROW NUMBER : ' 

MAXIMUM ROW NUMBER : ' 

MINIMUM COLUMN NUMBER : ' 

MAXIMUM COLUMN NUMBER : ' 

THE ELEVATION DATA FILE NAME : ' 

THE IMAGE DATA FILE NAME : ' 

THE I FRAME NUMBER OF IMAGE : ' 

THE J FRAME NUMBER OF IMAGE : ' 

WAIT APPROX. 1 MINUTE FOR SETUP*****' 



SUBROUTINE TER_CROP( IROW, LROW, ICOL, LCOL, IELEV2 ) 



THIS SUBROUTINE READS THE ORIGINAL TERRAIN GRID 
AND CONSTRUCTS A SMALLER GRID OF THE ELEVATION 
POINTS DESIRED. 



CHARACTER SWLAT* 8 , SWLON* 8 , DELLAT* 4 , DELLON* 4 , 
CHARACTER C0LS*4, ROWS*4 

INTEGER IELEV1( 210,239) , IELEV2( 50,50) , IROW, LROW 
INTEGER ICOL, LCOL, JV( 239) 

THE VALUES OF JV, SWLON, SWLAT DELLON, DELLAT, COLS 
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C AND ROWS ARE IGNORED. 

READ( 1,5) SWLON , SWLAT , DELLON , DELLAT , COLS , ROWS 
5 FORMAT( 1X,2( A8,2X) ,4( A4,2X) ) 

DO 20 M=1 ,238 

READ( 1,10) JV( M) , ( IELEV1( N, M) , N=1 , 210 ) 

10 FORMAT( I6,2x,20I6/( 8X,20I6) ) 

20 CONTINUE 

DO 40 N=IROW, LROW 
DO 30 M=ICOL, LCOL 

IELEV2( M+ 1- I COL, N+ 1- IROW ) = IELEV1( N, M ) 

30 CONTINUE 

40 CONTINUE 

RETURN 
END 
C 

Q ★★★★*★★★★★****★★★★**★****★★★★***★★★★★*★**★★★★**★*★★**★★★★★ 

SUBROUTINE READ I MAGE ( IMAGE ) 

C 

C THIS SUBROUTINE READS THE IMAGE DATA INTO AN ARRAY. 

C 

BYTE IMAGE( 512,512) 

C 

DO 10 IR=1 ,512 

READ( 4, REC=IR) ( I MAGE ( IC, IR) , IC=1,512) 

10 CONTINUE 

RETURN 
END 
C 

0 'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k’k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

SUBROUTINE TER_INTRP( IENDN, IENDM, IELEV2 ) 

C 

C THIS SUBROUTINE DETERMINES THE AVERAGE VALUE OF THE 

C ELEVATION DATA AND ASSIGNS THAT VALUE TO ANY UNKNOWN 

C POINTS 

C 

INTEGER IENDN, IENDM, IELEV2( 50, 50) 

DATA NEXT, I EL, IAVG/0,0,0/ 

C 

DO 20 N=l, IENDN 
DO 10 M=l, IENDM 
IF( IELEV2( M, N) . EQ. -32767)THEN 
GOTO 10 
ELSE 

NEXT=NEXT+ 1 
IEL=IEL+ IELEV2( M,N) 

END IF 

10 CONTINUE 

20 CONTINUE 

IAVG=IEL/NEXT 

C CHANGE UNKNOWN ELEVATION VALUES TO THE 

C CALCULATED AVERAGE. 

DO 40 N=l, IENDN 
DO 30 M=l, IENDM 
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IF( IELEV2( M,N) . EQ. -32767)THEN 
IELEV2( M, N ) =1 AVG 
END IF 

30 CONTINUE 

40 CONTINUE 
RETURN 
END 
C 

C ********************************************************* 
SUBROUTINE REAL_EL( IENDN, IENDM, IELEV2 ) 

C 

C THIS SUBROUTINE CREATS A REAL FILE OF THE ELEVATIONS 

C 

INTEGER IELEV2( 50,50) , IENDN, IENDM 
REAL AELEV( 239 ) 

C 

DO 20 N=l, IENDN 
K=0 

DO 10 M=1 , IENDM 
K=K+ 1 

AELEV( K) = IELEV2( M,N) 

10 CONTINUE 

WRITE( 2 , REC=N) ( AELEV( INT) , INT=1, IENDM) 

20 CONTINUE 

RETURN 
END 
C 

Q ********************************************************** 

SUBROUTINE TER_DMS( I ROW, LROW, I COL, LCOL, IENDM) 

C 

C THIS SUBROUTINE CONVERTS EACH ELEVATION POINT TO ITS 

C DMS EQUIVALENT. IT USES THE FACT THAT EACH ELEVATION 

C POINT REPRESENTS A ONE SECOND CHANGE IN LATITUDE 

C OR LONGITUDE FROM THE NEXT POINT. 

C REFERENCE POINT: LAT. 03 7 DEG. : 22 MIN. : 47 SEC. NORTH 

C LON. -122 DEG. : 05 MIN. : 03 SEC. WEST 

C OUR ELEVATION POINTS STAY WITHIN THE BOUNDS OF 037 

C DEGREES NORTH AND -122 DEGREES WEST, SO THESE VALUES 

C WILL EE ASSIGNED. 

C 

IMPLICIT DOUBLE PRECISION (A-Z) 

INTEGER IROW, LROW, ICOL, LCOL, IENDM, IL, JL 
REAL X , Y , Z , LATD , LATM , LATS , AELE V( 239) 

REAL LOND,LONM, LONS, HEIGHT 

PARAMETER( RLATM=22. ,RLATS=47. ,RLONM=5. ,RLONS=3. ) 

C 

LATD=37. 

LOND=- 122 . 

DO -20 I L= I ROW, LROW 
K=IL/60 
LATM=RLATM+ K 
LATS=RLATS + ( IL-K*60) 

IF( LATS. GE. 60. 0 ) THEN 
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LATS=LATS-60. 0 
LATM=LATM+ 1 . 0 
END IF 

I l=IL-( IROW-1) 

READ( 2 , REC=I 1 ) ( AELEV( I NT ) , I NT= 1 , I ENDM ) 

1=0 

DO 10 JL= I COL , LCOL 
K=JL/60 
LONM=RLONM-K 
LONS=RLONS-( JL-K*60) 

IF( LONS. LT. 0. 0 ) THEN 
LONM=LONM-l. 0 
LONS=LONS + 60. 0 
END IF 
LONM=-LONM 
LONS=-LONS 
1 = 1 + 1 

HEIGHT=AELEV( I) 

CALL DMS2XYZ( LAID, L ATM, L ATS, LOND, LONM, LONS, 



********************************************************** 



SUBROUTINE IM_REFIJ( IA, JA, I FRAME , JFRAME , I ENDM, IENDN) 

THIS SUBROUTINE CONSTRUCTS THE I AND J COORDINATE 
DATA FOR EACH ELEVATION POINT AND THEN CONVERTS THEM 
TO SCREEN COORDINATES AND STORES THEM IN ARRAYS 
I A AND JA. 

IMPLICIT DOUBLE PRECISION (A-Z) 

REAL OMEGA, PHI , KAPPA, X, Y, Z, XI , Y1 , Z1 , XO, YO 
REAL XIMA, YIMA A1 , A2 , B1 , B2 , Cl , C2 , FOCUS 
INTEGER I, J, I ENDM, IENDN, IA( 2500) ,JA( 2500) 

INTEGER IFRAME, JFRAME, I TOT 

DATA OMEGA, PH I, KAPPA/. 8341764,-. 4563699,3. 0761254/ 

DATA X0,Y0/0. 000002,0. 0/ 

DATA XI , Y1 , Z1/-2693765. 9,-4304520. 4,3859018. 3/ 

DATA Al, A2/20. 11323959,-6. 022849824/ 

DATA B1 , B2/6. 016207940,20. 10938801/ 

DATA Cl , C2/-34954. 59484, -22566. 71593/ 

DATA FOCUS/O. 153197/ 



HEIGHT, X, Y,Z) 



99 

10 

20 



WRITE( 3 , 99 )X, Y, Z 
FORMAT( 3( IX, F17. 7) ) 
CONTINUE 



CONTINUE 
ENDFILE( 3) 
REWIND( 3 ) 
RETURN 
END 



C 



ITOT=IENDN* IENDM 
DO 10 IR=1 , ITOT 
READ( 3,5, END=20 )X, Y, Z 
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FORMAT ( 3( IX, F17. 7) ) 

CALL PROJECT( X, Y, Z , OMEGA, PHI , KAPPA, XO , YO , XI , Y1 , Z1 , 
XIMA, YIMA, FOCUS) 

CALL XY2IJ( XIMA, YIMA, I , J, Al , A2 , B1 , B2 , Cl, C2 ) 

IA( IR)=I- I FRAME 
JA( IR) =4999- JFRAME- J 
10 CONTINUE 

20 REWIND( 3 ) 

RETURN 
END 

*1c*Jc1cit’k'k'k'k-k'k'k'k , k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'kjc'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'kje'k'k'k'k 

SUBROUTINE DMS2XYZ( LATD, LATM, LATS, LOND, LONM, LONS, 

LONS, HEIGHT, X,Y,Z) 

THIS SUBROUTINE CONVERTS DMS DATA TO X,Y, AND Z 
GEORECTANGULAR COORDINATES. 

IMPLICIT DOUBLE PRECISION (A-Z) 

REAL PHI , LAMDA , N , X , Y , Z, LATD, LATM, LATS 
REAL LOND, LONM, LONS, HEIGHT 
PARAMETER( PI =3. 14159265358793238) 

PARAMETER( Cl=180. ,C2=60. ,C3=3600. ) 

PARAMETER( E_SQUARE=0. 006768658 , A=6378206. 4) 

RAD I AN=P I /C 1 

PHI = ( LATD + LATM/C2 + L AT S/C 3 ) * RAD I AN 
LAMDA=( LOND+LONM/C2+LONS/C3 ) *RADIAN 
N=A/SQRT( 1-E_SQUARE*SIN( PHI )*SIN( PHI ) ) 

X=( N+HE IGHT ) *COS( PHI ) *COS( LAMDA) 

Y=( N+ HEIGHT ) *COS( PHI ) *SIN( LAMDA) 

Z=(N*( 1-E_SQUARE ) +HEIGHT ) *SIN( PHI) 

RETURN 
END 

********************************************************* 
SUBROUTINE PRO JECT( X , Y , Z , OMEGA , PHI , KAPPA , XO , YO , XI , 

Yl, Zl, XIMA, YIMA, FOCUS) 

THIS SUBROUTINE CONVERTS THE X,Y,Z GEORECTANGULAR 
COORDINATES TO XIMA AND YIMA WHICH ARE IMAGE 
PLANE COORDINATES. 

IMPLICIT DOUBLE PRECISION (A-Z) 

REAL Mil , M12 , M13 , M21 , M22 , M23 , M3 1 , M32 , M33 , DENOM 
REAL XIMA, YIMA, X, Y, Z, OMEGA, PHI , KAPPA 
REAL X0,Y0,X1,Y1,Z1, FOCUS 

Mll=COS( PHI ) *C0S( KAPPA) 

M12=COS( OMEGA) *SIN( KAPPA) + 

SIN( OMEGA) *SIN( PHI )*C0S( KAPPA) 

M13=SIN( OMEGA) *SIN( KAPPA)- 

COS( OMEGA) *SIN( PHI )*COS( KAPPA) 
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M2 1=-C0S( PHI )*SIN( KAPPA) 

M22=COS( OMEGA) *COS( KAPPA) - 

SIM( OMEGA) *SIN( PHI) *SIN( KAPPA) 

M23=SIN( OMEGA) *COS( KAPPA) + 

COS( OMEGA) *SIN( PHI ) *SIN( KAPPA) 

M31=SIN( PHI ) 

M32=-SIN( OMEGA) *COS( PHI ) 

M33=COS( OMEGA) *COS( PHI ) 

DENOM=M31*( X-Xl ) +M32*( Y-Yl ) +M33*( Z-Zl ) 

XIMA=XO- FOCUS* ( Mll*( X-Xl ) +M12*( Y-Yl ) +M13*( Z-Zl ) ) /DENOM 

YIMA=YO-FOCUS*( M21*( X-Xl ) +M22*( Y-Yl ) +M23*( Z-Zl ) ) /DENOM 

XIMA=XIMA* 1000000 

YIMA=YIMA* 1000000 

RETURN 

END 

********************************************************* 
SUBROUTINE XY2 1 J( XIMA, YIMA, I , J, Al, A2, Bl, B2 ,C1,C2 ) 

THIS SUBROUTINE TAKES THE IMAGE POINTS XIMA, YIMA 
AND CONVERTS THEM TO I,J ORIGINAL IMAGE COORDINATES. 

IMPLICIT DOUBLE PRECISION (A-Z) 

REAL XIMA, YIMA, A1 , A2 , B1 , B2 , Cl , C2 , DENOM 
INTEGER I,J 

DEN0M=A1*B2-B1*A2 

I=( ( XIMA- Cl ) *B2-( YIMA-C2)*B1) /DENOM 
J=-( ( XIMA- Cl ) *A2-( YIMA-C2 ) *A1 ) /DENOM 
RETURN 
END 

********************************************************* 
SUBROUTINE IM_REFAVG( IA, JA, IMAGE, IENDM, IENDN) 

THIS SUBROUTINE CONSTRUCTS THE FILE THAT IDENTIFIES 
THE ELEVATION POINTS THAT MAKE UP A TRIANGULAR PANEL 
AND ITS ASSOCIATED SAMPLED GREY SCALE VALUE. 

IMPLICIT DOUBLE PRECISION (A-Z) 

REAL SLOPE, YINT 
BYTE IMAGE( 512, 512) 

INTEGER I A( 2500) ,JA( 2500) , IY,M,N, IGREY1 , IGREY2 , L, LI 
INTEGER IENDM, I ENDN , NODE_A , NODE_B , NODE_C , NODE_D , IR 
INTEGER I COUNT 1, IC0UNT2, IT0T1, IT0T2 

DO 90 IR=1 , IENDN- 1 
IL=( IR-1)* IENDM 
DO 80 N=1+IL, IENDM-l+IL 
NODE_A=N 
NODE_B=N+ 1 
NODE C=N+ IENDM 
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SLOPE=( JA( NODE_C) - JA( NODE_B) ) *1. 0/ 

( IA(NODE_C)-IA(NODE_B) )*1. 0 
YINT=1. 0* JA( N0DE_B ) -SLOPE* IA( N0DE_B ) 

IT0T1=0 
IT0T2=0 
I COUNT 1=0 
IC0UNT2=0 

DO 70 M=IA( N0DE_B) , IA( N0DE_A) 

IY=( SLOPE*M+YINT) 

DO 50 L=JA( NODE_B) , IY 
IT0T1=IT0T1+ IMAGE( M, L) 

I COUNT 1= I COUNT 1+1 
50 CONTINUE 

IF( M. LT. IA( N0DE_A) )THEN 
DO 60 L=IY+1, JA(N0DE_C) 

I TOT2 = I T0T2 + I MAGE ( M , L ) 

ICOUNT2=ICOUNT2+ 1 
60 CONTINUE 

END IF 

70 CONTINUE 

Ll=( N-( IR-1) ) *2 
IGREY1=IT0T1/IC0UNT1 
IGREY2=ITOT2/ICOUNT2 
NODE_D=NODE_C+ 1 

WRITE( 20,REC=L1-1)NODE_A,NODE_B,NODE_C, IGREY1 
WRITE( 20, REC=L1 )NODE_B, N0DE_D , NODE_C , IGREY2 
80 CONTINUE 

90 CONTINUE 
RETURN 
END 

-k-k’k'k'k'k'k'k-k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k-kic'k-k-k-k’k’k’k’k'k’k’k-k-k'k-k-k'k’k’k'k-k-k-k'k'k'k 

SUBROUTINE EXT_0RIEN( Mil , M12 , M13 , M21 , M22 , M2 3 , 

M31 , M32 , M33 , IENDM, IENDN) 

THIS ROUTINE DETERMINES THE M MATRIX PARAMETERS 
FOR ROTATION OF THE IMAGE PLANE TO DESIRED 
LOCATION FOR VIEWING IN THE SYNTHESIZED IMAGE. 

IMPLICIT DOUBLE PRECISION (A-Z) 

REAL MAGN_X , MAGN_Y , MAGN_Z , X , Y , Z 
REAL X_CORD( 2500 ) , Y_CORD( 2500) , Z_CORD( 2500) 

REAL X_VECX, X_VECY, X_VECZ , Y_VECX, Y_VECY, Y_VECZ 
REAL Z_VECX,Z_VECY,Z_VECZ 
INTEGER IENDM, IENDN, ITOT, IP, IR 
C 

I TOT= I ENDM* I ENDN 
DO 10 IR=1 , ITOT 
READ( 3,5, END=20 )X, Y, Z 
X_C0RD( IR ) =X 
Y_CORD( IR)=Y 
Z_CORD( I R ) =Z 

5 FORMAT( 3( IX, FI 7. 7 ) ) 
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10 CONTINUE 

20 REWIND( 3 ) 

Y_VECX=-( X_C0RD( 1 ) ) 

Y_VECY=-( Y_C0RD( 1 ) ) 

Y_VECZ=-( Z_C0RD( 1 ) ) 

I P= I TOT- I ENDM+ 1 
Z_VECX=X_CORD( 1 ) -X_C0RD( IP ) 

Z_VECY=Y_CORD( 1 ) -Y_C0RD( IP ) 

Z_VECZ=Z_C0RD( 1 ) -Z_C0RD( IP) 

USE THE CROSS PRODUCT OF Y CROSS Z TO OBTAIN 
THE X VECTOR. 

X_VECX=( ( Y_VECY*Z_VECZ ) -( Y_VECZ*Z_VECY) ) 

X_VECY=( ( Y_VECZ*Z_VECX)-( Y_VECX*Z_VECZ ) ) 

X_VECZ=( ( Y_VECX* Z_VECY ) - ( Y_VECY*Z_VECX ) ) 

MAGN_Z=SQRT( ( Z_VECX**2 ) +( Z_VECY**2 ) +( Z_VECZ**2 ) ) 
MAGN_X=SQRT( ( X_VECX**2 )+( X_VECY**2 ) +( X_VECZ**2 ) ) 
MAGN_Y=SQRT( ( Y_VECX**2 ) + ( Y_VECY**2 ) + ( Y_VECZ**2 ) ) 

Mil =X_VECX/MAGN_X 
M 1 2 =X_ VEC Y/MAGN_X 
Ml 3 =X_VECZ/MAGN_X 
M2 1 =Y_VECX/MAGN_Y 
M22=Y_VECY/MAGN_Y 
M2 3 =Y_VECZ/MAGN_Y 
M3 1=Z_VECX/MAGN_Z 
M3 2 = Z_VEC Y/MAGN_Z 
M3 3 =Z_VECZ/MAGN_Z 
RETURN 
END 

SUBROUTINE AFFIN( A1 , A2 , B1 , B2 , Cl , C2 ) 

THIS SUBROUTINE ASSIGNS OR CALCULATES THE 
COEFFICIENTS TO BE UTILIZED IN THE TRANSFORM FROM 
IMAGE COORDINATES TO SCREEN COORDINATES. 

IMPLICIT DOUBLE PRECI SI0N( A-Z ) 

REAL A1,A2,B1,B2,C1,C2,XIMA_MAX,YIMA_MAX, I_MAX, J_MAX 
DATA I_MAX, J_MAX/512. 0,512. 0/ 

XIMA_MAX=1600. 0 
YIMA_MAX=1600. 0 
C 1 =X I MA_MAX 
C2=0. 0 
A2=0. 0 
B1=0. 0 

Al=-XIMA_MAX/( I_MAX*1. 0) 

B2=YIMA_MAX/( J_MAX*1. 0) 

RETURN 
END 

********************************************************** 
SUBROUTINE OBS_LOC( XI , Y1 , Z1 , FOCUS ) 
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c 

C THIS SUBROUTINE CALCULATES THE NEW OBSERVER X1,Y1,Z1 

C LOCATION FROM DESIRED LAT. AND LONG. INPUTS AS WELL 

C AS PROVIDE THE FOCAL LENGTH. 

C 

IMPLICIT DOUBLE PRECISION (A-Z) 

REAL LATD , L ATM , LAT S , LOND , LONM , LONS , HE I GHT 
REAL X1 / Y1,Z1,X,Y / Z / FOCUS 
C 

LATD=37. 0 
WRITE( 6, *) ' INPUT 
READ( 5,5) LATM 
WRITE( 6 , * ) ' 

READ( 5,5) LATS 
LOND=-122. 0 
WRI TE( 6, * ) ' INPUT 
READ( 5,5) LONM 
WR I TE ( 6, *) ' 

READ( 5, 5) LONS 
5 FORMAT( F5 . 1 ) 

WRITE( 6,*) ' INPUT 
READ( 5, 10) HE I GHT 
10 FORMAT( F6. 1 ) 

FCCUS=0. 015 

CALL DMS2XYZ( LATD, LATM, LAT S , LOND , LONM , LONS, HEIGHT, 

X, Y, Z ) 

X1=X 

Y1=Y 

Z1=Z 

RETURN 

END 

C 

Q ^^^^^^^^'k’k’k'k'k-k-k’k’k'k’k’k'k-k'k’k’k'k-k'k'k'k’k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k’k'k'k'k'k 

SUBROUTINE NEW_I J( XI , Y1 , Z1 , FOCUS , A1 , A2 , B1 , B2 , Cl , C2 , 

IENDM, IENDN, IA, JA) 

C 

C THIS SUBROUTINE COMPUTES THE NEW IA AND JA SCREEN 

C COORDINATES FROM THE GIVEN OBSERVER LOCATION. 

C 

IMPLICIT DOUBLE PRECISION (A-Z) 

REAL XI , Y1 , Z1 , FOCUS , Mil , M12 , M13 , M2 1 , M22 , M23 , M3 1 
REAL XO , YO , X, Y, Z, XIMA, YIMA, A1 , A2 , B1 , B2 , Cl , C2 
REAL M3 2 , M3 3 , DENOM 

INTEGER IENDM, IENDN, ITOT, IR, IA( 2500 ) , JA( 2500 ) , I , J 
DATA XO , Y0/0. 0,0. 0/ 

C 

I TOT= IENDN* IENDM 
DO 10 IR=1 , 2500 
I A( IR)=0 
JA( IR) =0 
10 CONTINUE 

CALL EXT_ORIEN( Ml 1 , M12 , M13 , M2 1 , M22 , M23 , M3 1 , M32 , M33 , 

IENDM, IENDN) 



OBSERVER LATITUDE IN-MINUTES( REAL ) : ' 

-SECONDS (REAL): ' 

OBSERVER LONGITUDE IN-MINUTES( REAL ) : ' 

-SECONDS( REAL): ' 

OBSERVER HEIGHT-METERS( REAL) : ' 
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DO 20 IR=1 , I TOT 
READ( 3, 15,END=30)X,Y,Z 
15 FORMAT( 3( IX, F17. 7 ) ) 

DENOM=M31*( X-Xl )+M32*( Y-Yl ) +M33*( Z-Zl ) 

XIMA=XO-FOCUS*( M11*(X-X1)+M12*( Y-Y1)+M13*( Z-Zl) )/ 
DENCM 

Yl MA=Y0- FOCUS *( M21*( X-Xl ) +M22*( Y-Yl ) +M23*( Z-Zl ) )/ 
DENOM 

XIMA=X IMA* 1000000. 0 
YIMA=YIMA* 100 0000. 0 

CALL XY2 1 J( XIMA, YIMA, I , J, A1 , A2 , B1 , B2 , Cl , C2 ) 

IA( IR)=I 
JA( IR)=J 
20 CONTINUE 

30 REWIND( 3 ) 

RETURN 

END 

C 

Q : k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k’k'k'k'k'k'k'k'k'k-k'k'k'k'k'k 

SUBROUTINE N0DE_DPTH( XI , Yl , Z1 , DEPTH, IENDM, IENDN) 

C 

C THIS SUBROUTINE CALCULATES THE DISTANCE OF THE 

C TRANSLATED ELEVATION POINTS TO THE NEW OBSERVER 

C LOCATION. 

C 

IMPLICIT DOUBLE PRECISION A-Z ) 

REAL XI, Yl, Z1,DEPTH( 2500) 

INTEGER IENDM, IENDN, IR, ITOT 
C 

I TOT= IENDM* IENDN 
DO 10 IR=1 , ITOT 
READ( 3,5, END=20 )X, Y, Z 
5 FORMAT( 3( 1X,F17. 7) ) 

DEPTH( IR)=SQRT( ( ( Xl-X) **2 ) + ( ( Yl-Y)**2)+( ( Zl-Z)**2) ) 

10 CONTINUE 

20 REWIND( 3 ) 

RETURN 

END 

C 

C ****** **************************************************** 

SUBROUTINE FILL( IA, JA, DEPTH, IMAGE, IENDM, IENDN) 

C 

C THIS SUBROUTINE DETERMINES THE HIDDEN SURFACES AND 

C CONSTRUCTS THE TRANSLATED IMAGE. A Z_BUFFER IS USED 

C TO HOLD THE COMPUTED DEPTHS FROM THE OBSERVER TO THE 

C GEOGRAPHIC POSITION. A PLANE EQUATION IS CONSTRUCTED 

C FROM THREE ELEV. POINTS. THIS EQUATION IS USED TO 

C DETERMINE THE DEPTH OF ALL POINTS WITHIN THE PLANE. 

C 

IMPLICIT DOUBLE PRECISION A-Z) 

BYTE I MAGE ( 512,512) 

REAL DEPTH( 2500) , Z_BUFF( 512 , 512 ) , XI , X2 , X3 , Yl , Y2 , Y3 
REAL Z1 , Z2 , Z3 , Z_DPTH, Cll , C12 , C13 , C21 , C22 , C23 , C31 
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REAL C32 , C33 , DET , A_COEF , B_COEF , C_COEF 
INTEGER I GREY, NODE_A, NODE_B , NODE_C, I_MIN, I_MAX, J_MIN 
INTEGER J_MAX, FRAME( 512,512 ) , IA( 2500) , JA( 2500) , IENDM 
INTEGER IENDN, IR, I, J,K, L, IPLANES 
C 

C FIRST DETERMINE THE NUMBER OF PLANES AND INITIALIZE 

C THE Z_BUFFER AND IMAGE TO 0. 

C 

IPLANES=( ( IENDM-1 ) *2 ) *( IENDN- 1 ) 

DO 20 K=1 , 512 
DO 10 L=1 , 512 
I MAGE ( L, K)=0 
Z_BUFF( L,K)=0 
10 CONTINUE 

20 CONTINUE 

C 

C DETERMINE THE COEFFICIENTS OF THE PLANE EQUATION 

C FROM THE THREE ELEVATION POINTS. 

C 

DO 70 IR=1, IPLANES 

READ( 2 0 , REC= I R ) NODE_A , N0DE_B , N0DE_C , I GREY 
Xl=IA(NODE_A) 

Yl= JA( NODE_A ) 

Z1=DEPTH( NODE_A) 

X2=IA(NODE_B) 

Y2=JA( NODE_B ) 

Z2=DEPTH(NODE_B) 

X3=IA( NODE_C) 

Y3=JA( NODE_C ) 

' Z3=DEPTH( NODE_C) 

C DETERMINE THE COFACTOR ELEMENTS 

Cll=( ( Y2*Z3)-( Y3*Z2) ) 

C12=-( ( X2*Z3 ) -( X3*Z2 ) ) 

C13=( ( X2*Y3 ) -( X3*Y2 ) ) 

C21=-( ( Yl*Z3)-( Y3*Z1) ) 

C22=( (X1*Z3)-(X3*Z1) ) 

C23=-( (X1*Y3)-(X3*Y1) ) 

C31=( ( Yl*Z2)-( Y2*Z1) ) 

C32=-( (X1*Z2)-(X2*Z1) ) 

C33 = ( (X1*Y2)-(X2*Y1) ) 

C CALCULATE THE DETERMINANT 

C 

DET=(Xl*Cll)+( Yl*C12)+( Z1*C13) 

C 

C THE COEFFICIENTS ARE DETERMINED FROM MULTIPLYING 

C THE reciprocal of the determenant with the 

C transpose of the cofactors called the adjoint. 

C 

A_C0EF=-( ( C11+C21+C31 )/DET) 

B_C0EF=-( ( C12+C22+C32 )/DET) 

C_COEF=-( ( C13+C23+C33 )/DET ) 

IF( C_C0EF. EQ. 0. 0)GOTO 70 

C DETERMINE THE MAXIMUM AND MINIMUM VALUES TO TEST 
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C FOR THE FILL ALGORITHUM. 

C 

I_MAX=MAX( I A( NODE_A ) , IA( NODE_B ) , IA( NODE_C ) ) 

I_MIN=MIN( IA(NODE_A) , IA(NODE_B) , IA(NODE_C) ) 

J_MAX=MAX( JA( NODE_A ) , JA( NODE_B ) , JA( NODE_C ) ) 

J_MIN=MIN( JA( NODE_A ) , JA( NODE_B ) , JA( NODE_C ) ) 

I F( I_MIN. LT. 1. OR. I_MAX. GT. 512 )GOTO 70 
I F( J_MIN. LT. 1. OR. J_MAX. GT. 512) GOTO 70 
C 

C CLEAR THE REFERENCE FRAME BUFFER AND CALL THE 

C FRAME FILL SUBROUTINE. 

C 

DO 40 L=I_MIN, I_MAX 
DO 30 K=J_MIN, J_MAX 
FRAME ( L / K)=0 
30 CONTINUE 

40 CONTINUE 

CALL FRAME_F I L( NODE_A , NODE_B , NODE_C , FRAME , IA, JA, 

J_MAX , J_M I N ) 

DO 60 J=J_MIN, J_MAX 
DO 50 I=I_MIN, I_MAX 
I F( FRAME( I , J ) . EQ. 1 ) THEN 

Z_DPTH=-( l + ( A_COEF*I ) + ( B_COEF*J) )/C_COEF 
IF( Z_BUFF( I , J) . EQ. 0. 0. OR. Z_DPTH. LT. Z_BUFF( I , J) )THEN 
Z_BUFF( I, J)=Z_DPTH 
I MAGE ( I , J)=IGREY 
END IF 
END IF 

50 CONTINUE 

60 CONTINUE 

70 CONTINUE 

DO 90 J=1 , 512 

WRITE( 2 1 , REC=J ) ( I MAGE ( I , J ) , 1=1 , 512 ) 

90 CONTINUE 

RETURN 
END 
C 

Q 'k'k’k’k'k'k'k'k'k'k'k’k'k'k'k'k'k'k'k'k'k'k'k'k'k'k’k'k’k'k’k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 

SUBROUTINE FRAME_FIL( NODE_A, NODE_B , NODE_C , FRAME , 

I A , JA , J_MAX , J_M I N ) 

C 

C THIS SUBROUTINE CONSTRUCTS AN EDGE LIST FROM THE 

C THREE NODES PASSED IN THE ROUTINE. THESE EDGES 

C ARE USED IN A POLYGON FILL ROUTINE USEING A FRAME 

C BUFFER AND Y_BUCKET 

IMPLICIT DOUBLE PRECI SION( A-Z ) 

REAL I_INCPT, DELTA_I , DELTA_J, AEL( 3,4) 

REAL XN0DE1 , XN0DE2 , DX 

INTEGER NODE_A , NODE_B , NODE_C , FRAME ( 512,512), IA( 2500) 
INTEGER JA( 2500 ) , J_MAX, J_MIN, NODE( 4) , I S , NODE 1 , N0DE2 
INTEGER N_H I GH , N_LOW , IT, IU, J_BUCKET( 512) , IR, II, 12 
INTEGER ICNT, LCNT 
C 
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DO 10 IS=1 , 512 
J_BUCKET( IS)=0 
10 CONTINUE 

DO 30 IS=1 , 3 
DO 20 IR=1 , 4 
AEL( IS, IR)=0. 0 
20 CONTINUE 

30 CONTINUE 

NODE( 1 ) =NODE_A 
NODE( 2 ) =NODE_B 
NODE( 3 ) =NODE_C 
NODE( 4 ) =NODE_A 
DO 40 I S=1 , 3 
NODEl=NODE( IS) 

NODE2=NODE( IS+1) 

IF( JA( NODE1 ) . GE. JA( NODE2 ) ) THEN 
N_H I GH=NODE 1 
N_LOW=NODE2 
ELSE 

N_HIGH=NODE2 
N_LOW=NODE 1 
END IF 

I_INCPT=IA( N_HIGH) 

DELTA_J=( JA( N_HIGH) - JA( N_LOW) ) 

IF( DELTA_J. EQ. 0. 0 ) THEN 
GOTO 40 
ELSE 

DELTA_I=-( IA(N_HIGH)-IA(N_LOW) )/DELTA_J 
END IF 

IT=JA(N_HIGH) 

IU=J_BUCKET( IT) 

IF( IU. EQ. 0 ) THEN 
J_BUCKET( IT) = IS 
ELSE 

AEL( IU, 4)=IS 
END IF 

AEL( IS, 1 )=I_INCPT 
AEL( I S , 2 ) =DELTA_I 
AEL( I S , 3 ) =DELTA_J 
40 CONTINUE 

IT=J_BUCKET( J_MAX) 

IU=INT( AEL( IT, 4) ) 

XN0DE1=AEL( IT, 1) 

XNODE2=AEL( IU, 1) 

DX=XN0DE1-XN0DE2 
IF(DX. LE. 0. 0 ) THEN 
1 1=NINT( XN0DE1 ) 

I2=NINT( XN0DE2 ) 

ELSE 

1 1=NINT( XN0DE2 ) 

I2=NINT( XN0DE1 ) 

END IF 

DO 50 I R= I 1 , 12 
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FRAME( I R , J_MAX ) = 1 
50 CONTINUE 

AEL( IT / 3)=AEL( 11,3)-!. 0 
AEL( IU, 3)=AEL( IU,3)-1. 0 
I CNT= J_MAX - 1 
LCNT=J_MIN 

DO 70 IS=ICNT / LCNT / -1 
I F( J_BUCKET( I S ) . EQ. 0 ) THEN 
XN0DE1=XN0DE1+AEL( IT, 2) 

XNODE2=XNODE2+AEL( IU,2) 

ELSE IF( AEL( IT, 3 ). LE. 0. 0)THEN 
I T= J_BUCKET( IS) 

XN0DE1=AEL( IT, 1) 

XNODE2=XNODE2+AEL( IU,2) 

ELSE 

IU=J_BUCKET( IS) 

XN0DE1=XN0DE1+AEL( IT, 2) 

XN0DE2=AEL( IU, 1) 

END IF 

DX=XNODE 1 -XN0DE2 
I F( DX. LE. 0. 0 ) THEN 
I 1 =N I NT ( XNODE 1 ) 

I2=NINT( XN0DE2 ) 

ELSE 

I 1=NINT( XN0DE2 ) 

I2=NINT( XNODE 1 ) 

END IF 

DO 60 IR=I 1 , 12 
FRAME ( IR, IS )=1 
60 CONTINUE 

AEL( IT, 3 ) =AEL( IT,3)-1. 0 
AEL( IU, 3 ) =AEL( IU,3)-1. 0 
70 CONTINUE 
RETURN 
END 
C 

0 'k'k'k'k'k'k'k'k'k’k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k’k'k'k'k'k'k’k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k'k 
0 Tfr******************************************************** 
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